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Abstract 

This  thesis  provides  a  theoretical  framework  for  haptics,  the  study  of  exploration  and  ma¬ 
nipulation  using  hands.  Be  it  human  or  robotic  research,  an  understanding  of  the  nature  of 
contact,  grasp,  exploration,  and  manipulation  is  of  singular  importance.  In  human  haptics  the 
objective  is  to  understand  the  mechanics  of  hand  actions,  sensory  information  processing,  and 
motor  control.  While  robots  have  lagged  behind  their  human  counterparts  in  dexterity,  recent 
developments  in  tactile  sensing  technology  have  made  it  possible  to  build  sensor  arrays  that  in 
some  way  mimic  human  performance.  We  believe  that  a  computational  theory  of  haptics,  that 
investigates  what  kind  of  sensory  information  is  necessary  and  how  it  has  to  be  processed  is 
beneficial  to  both  human  and  robotic  research. 

Human  and  robot  tactile  sensing  can  be  accomplished  by  arrays  of  mechanosensors  embed¬ 
ded  in  a  deformable  medium.  When  an  object  comes  in  contact  with  the  surface  of  the  medium, 
information  about  the  shape  of  the  surface  of  the  medium  and  the  force  distribution  on  the 
surface  is  encoded  in  the  sensor  signals.  The  problem  for  the  central  processor  is  to  reliably 
and  efficiently  infer  the  object  properties  and  the  contact  state  from  these  signals.  In  the  first 
part  of  the  thesis  we  discuss  the  surface  signal  identification  problem:  the  processing  of  sensor 
signals  resulting  in  algorithms  and  guidelines  for  sensor  design  that  give  optimal  estimates  of 
the  loading  and  displacement  distributions  on  the  surface  of  the  fingerpad. 

In  the  second  part  of  the  thesis  we  focus  on  how  the  information  obtained  from  such  optimal 
sensing  can  be  used  for  exploration  of  objects.  We  argue  that  an  accurate  reconstruction 
of  object  properties  can  occur  using  two  basic  building  blocks  of  Exploration  Strategy  and 
Finger  Control.  Exploration  Strategy  pertains  to  the  problem  of  inferring  object  properties 
such  as  shape,  texture  and  compliance  and  interface  properties  such  as  state  of  contact,  from 
the  estimated  surface  signals.  This  involves  determining,  in  each  case,  what  kind  of  sensor 
information  and  what  kind  of  action  is  needed.  Finger  Control  refers  to  the  transformation  of  the 
action  needed  into  a  command  trajectory  for  the  fingerpad,  which  defines  the  desired  direction 
of  movement  for  manipulation.  We  define  and  anlyze  the  components  of  both  these  blocks, 
provide  explicit  mathematical  formulation,  and  give  numerical  examples  where  appropriate. 
Our  formulation  of  this  computational  theory  of  haptics  is  independent  of  implementation  so 
that  it  is  applicable  to  both  robots  and  humans. 

Thesis  Supervisor:  Anuradha  M.  Annaswamy 
Title:  Associate  Professor 
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Chapter  1 


Introduction 


Compared  to  vision  and  hearing,  we  might  attach  less  importance  to  the  sense  of  touch,  pri¬ 
marily  because  the  former  senses  can  easily  be  temporarily  occluded  and  we  are  therefore  more 
aware  of  their  presence.  We  therefore  tend  to  underestimate  the  importance  of  the  sense  of 
touch  but  on  second  thought  we  realize  that  without  it,  our  life  would  be  much  harder  and  pos¬ 
sibly  more  difficult  than  without  vision  or  hearing  [48].  We  apply  touch  for  almost  any  activity 
involving  direct  interaction  with  our  environment  and  in  robotic  research  it  has  been  observed 
that  controlling  this  interaction  is  very  difficult  without  tactile  information.  The  research  into 
the  nature  of  tactile  sensing  has  been  conducted  more  or  less  separately  in  the  fields  of  robotics 
and  physiology,  but  the  primary  difference  between  these  fields  lies  only  the  implementation. 
We  however  believe  that  a  common  framework  and  theory  can  be  constructed,  revealing  the 
limitations  and  possibilities  of  touch  and  manipulation  and  further  explicitly  formulating  the 
information  processing  involved. 

The  use  of  the  word  Haptics  used  in  the  title  of  this  thesis  is  primarily  confined  to  a  circle 
of  people  working  on  tactile  sensing  research  and  in  other  related  areas.  The  term  refers  to 
the  mechanics  and  the  information  processing  involved  in  tactual  exploration,  manipulation 
and  perception,  and  it  is  used  in  robotic  as  well  as  physiological  circles.  A  computational 
theory  of  haptics  is  then  a  formulation  of  the  relevant  information  transfer  and  processing  and 
will  necessary  include  mechanistic  as  well  as  signal  processing  models.  This  formulation  is  the 
primary  subject  of  this  thesis.  The  problems  we  concentrate  on,  are  (i)  the  mechanical  modeling 
of  compliant  fingerpads,  (ii)  the  information  processing  of  signals  from  strain  or  stress  sensors 
embedded  in  the  fingerpad,  (iii)  the  inference  of  object  and  environment  properties  from  these 
signals,  and  finally  (iv)  the  relation  of  the  above  to  the  movement  of  hands  and  fingers  during 
haptic  exploration. 

The  contributions  of  this  thesis  can  be  summarized  as  the  presentation  of  a  computational 
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theory  of  touch,  which  we  can  itemize  as: 


1.  A  complete  solution  for  a  3D  linear  model  of  the  fingerpad  with  general  normal  and  shear 
loading  distributions,  also  including  elastic,  viscous  and  inertial  effects. 

2.  Transfer  function  formulation  of  the  above  solution,  an  approach  that  provides  access  to 
well  developed  and  effective  signal  processing  tools  such  as  the  Fast  Fourier  Transform 
(FFT)  and  the  Kalman  and  Wiener  filters  which  we  evoke  to  decode  surface  signals  from 
subsurface  sensor  signals. 

3.  We  present  an  optimal  set  of  sensors  selected  from  all  linear  strains  and  stresses.  So  as 
to  obtain  maximum  spatial  bandwidth  in  the  sensor  signals  without  a  directional  bias. 

4.  A  detailed  framework  delineating  definitions  and  providing  consistent  formulation  of  ex¬ 
ploratory  procedures  and  control  strategies  that  can  be  used  for  the  identification  of  object 
properties  such  as  shape,  texture  and  compliance. 

The  thesis  is  composed  of  four,  largely  independent  chapters  in  addition  to  appendices, 
each  discussing  a  specific  aspect  of  the  above  problem.  The  chapters  are  designed  so  that 
they  can  read  as  independent  papers.  Therefore,  a  more  detailed  introduction  and  background 
information  such  as  literature  review,  are  to  be  found  in  the  introduction  part  of  each  chapter, 
and  we  refer  the  reader  to  those  sections  for  that  information.  Chapter  2  is  slightly  different, 
in  that  it  concentrates  on  giving  tutorial  information,  primarily  on  the  mathematical  methods 
and  tools  used  in  subsequent  chapters.  Some  parts  of  that  chapter  may  therefore  be  repeated 
later  in  the  thesis,  but  then  in  order  to  preserve  the  continuity  within  and  independence  of  each 
chapter. 

Chapterwise  the  thesis  is  organized  as  follows:  Chapter  2  contains  a  tutorial  on  the  math¬ 
ematical  methods  applied  in  the  subsequent  chapters.  In  chapter  3  we  present  a  model  valid 
under  static  or  quasi-static  conditions  and  where  we  further  show  how  our  methods  can  be 
applied  to  predict  the  onset  of  slip  before  the  finger  starts  to  slide.  In  chapter  4  we  delineate  a 
dynamic  version  of  the  fingerpad  model  including  viscoelastic  and  inertial  effects.  We  further 
analyze  the  solution  and  derive  conditions,  under  which  the  inertial  effects  can  be  ignored  but 
that  considerably  simplifies  the  solution.  In  chapter  5  we  discuss  exploration  strategies,  con¬ 
necting  object  properties,  sensors,  signal  processing  and  hand  or  fingerpad  movements.  Finally, 
in  appendices  A-F  we  give  detailed  derivations  as  well  as  some  mathematical  manipulations 
and  observations  regarding  the  solutions  presented  in  the  chapters  delineated  above. 
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Chapter  2 


Mathematical  Preliminaries 


1.  Introduction 

Most  people  working  on  mechanical  systems  are  aware  of  integral  transforms,  such  as  Fourier 
and  Laplace  transforms,  as  tools  used  to  solve  differential  equations  or  the  analysis  of  measure¬ 
ment  data.  The  general  attitude  however  differs  widely,  depending  mostly  on  the  discipline 
people  belong  to.  Control  engineers  work  extensively  with  integral  transforms  and  are  in  gen¬ 
eral  quite  comfortable  in  expressing  the  solutions  in  the  form  of  a  transform.  In  many  cases  the 
answer  would  be  a  transfer  function  describing  a  system  without  the  input  being  prescribed.  A 
mechanics  engineer,  on  the  other  hand,  is  more  likely  to  bring  the  solution  one  step  further  by 
prescribing  some  specific  input  (boundary  conditions)  and  then  go  through  effort  of  obtaining 
the  actual  solution  to  that  input.  Both  approaches  have  their  merit,  but  ideally  one  would 
apply  both  and  use  specific  solutions  as  well  as  transfer  functions  to  explain  or  predict  the 
systems  response  to  generic  inputs. 

In  this  chapter  we  will  attempt  to  bring  these  two  above-mentioned  disciplines  closer  to¬ 
gether  in  the  application  of  integral  transforms  to  contact  problems.  We  will  first  give  the  basic 
theory  of  integral  transforms  and  then  apply  it  to  several  structures  of  differing  complexity. 
The  structures  we  will  analyze,  will  vary  from  the  statically  loaded  spring  (Winkler)  founda¬ 
tion  to  the  dynamic  loading  of  the  semi-infinite  elastic  halfspace.  By  going  from  the  simple 
structures  to  the  more  complex  and  by  looking  at  static  loading  before  dynamic,  we  hope  to 
introduce  various  modeling  details  and  their  consequences  separately.  Through  the  analysis  of 
the  transfer  function  solution  for  each  structure  and  by  comparing  that  solution  to  some  of  the 
available  solutions  we  hope  that  in  the  end,  the  reader  will  be  well  acquainted  and  comfortable 
with  both  forms  of  solutions. 
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The  use  of  integral  transforms  in  contact  mechanics  has  mainly  been  promoted  by  Sneddon 
and  for  a  detailed  discussion  cf.  [70,  71],  but  more  applications,  models,  and  solution  methods 
are  given  in  [15,  61,  51].  Finally,  for  a  general  reference  on  contact  mechanics  cf.  [36],  but  on 
viscoelasticity  cf.  [83]. 

This  chapter  is  organized  as  follows:  In  section  2,  we  give  the  definition  and  some  basic 
properties  of  Fourier  (ID  and  2D),  and  Laplace  transforms.  We  also  discuss  briefly,  some 
variants  of  the  Fourier  transforms,  namely  the  finite  and  discrete  Fourier  transforms.  Section 
3  is  devoted  to  static  contact  problems,  and  section  4  to  dynamic  problems.  We  conclude  by 
summarizing  our  conclusions  in  section  5. 


2.  Integral  Transforms 


The  use  of  integral  transforms  for  the  solution  of  partial  differential  equations  (PDE)  is  based 
on  the  method  of  separation  of  variables  introduced  by  d’Alembert,  Bernoulli  and  Euler  in  the 
mid-eighteenth  century  [71].  The  transform  concept  becomes  convenient  when  the  solutions  are 
parameterized  by  the  separation  constants  that  characterize  the  method. 


Example:  When  using  the  method  of  separation  of  variables  to  obtain  a  solution  of  Laplace’s 
equation  in  Cartesian  coordinates 


d2<p  d2(f>  d24> 
dx 2  dy 2  +  dz2 


(2.1) 


the  usual  procedure  is  to  assume  that  the  solution  can  be  expressed  as 


4  =  X(x)Y(y)Z(z). 


Substituting  this  into  Eq.  (2.1)  and  dividing  by  X(x)Y(y)Z(z)  we  obtain 

X"(x)  ,  Y"(y)  ,  Z"{z) 


+ 


+ 


X(x)  Y(y)  Z{z) 


=  0. 


We  observe  that  the  first  term  is  a  function  of  x  only,  the  second  of  y  only  and  the  third  of  only 
z.  They  must  therefore  each  be  a  constant.  Let 

X"  Y"  Z" 

jr  =  &X  ~p~  ~  tty  then  ^  Cly , 

where  ax  and  ay  can  be  any  constants,  complex  or  real.  We  then  have  that 

Y  _  e±ia:\/S7  y  —  e±??/\/“7  2  — 
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and  that 


(j)  =  e  ±jxs/a^±jy^;  e±Zy/ax+ay 


is  a  solution  of  Eq.  (2.1)  for  all  ax,  ay. 

This  solution  method  can  be  applied  to  a  wide  variety  of  PDEs  and  we  will  now  discuss 
how  it  can  generally  be  developed  to  solve  PDEs  with  prescribed  boundary  conditions. 


Suppose  that  we  have  to  determine  a  displacement  field  u,  the  solution  of  a  homogeneous 
PDE  of  the  form 


L  [u(r)]  —  0  r  £  D 


in  which  r  =  (aq,  aq,  aq)  is  the  position  vector  of  a  field  point  in  a  domain  D,  and  L[-\  is  a 
linear  differential  operator  in  the  variables  aq,  aq,  and  23.  The  set  which  the  domain  D  is  a 
part  of  is  in  most  cases  prescribed  by  the  geometry  and  we  will  assume  here,  for  the  sake  of 
clarity,  that  D  C  H3  and  the  operator  therefore  given  in  Cartesian  coordinates.  The  discussion 
given  here  will  also  apply  to  other  coordinate  systems,  such  as  cylindrical,  spherical  and  others 
and  we  will,  where  applicable,  give  special  comments  regarding  to  them. 


Assuming  that  the  method  of  separation  of  variables  enables  us  to  find  a  solution  ??(?•,  aq ,  aq) 
involving  2  arbitrary  (separation)  constants  (aq,  aq),  whose  values  are  taken  from  a  set  Q.  Since 
the  operator  L[-\  is  linear,  the  function 

u(r)  —  /  ?;(r,  aq,  aq)F(aq,  aq^aqdaq 
J  n 


will  also  be  a  solution  for  a  class  of  functions  F(-).  If  we  in  addition  have  the  boundary 
condition1 

u(.zq,  ifq,  23  A3)  /(x^,  aq ) 

A'3  being  a  constant,  then  we  must  choose  f(a q,aq)  such  that 

/(aq,aq)  =  /  A’ (aq,  aq,  aq,aq)/(aq,  a2)daqdaq  (aq,  aq)  £  D  fl  {aq  =  A"3} 
in 


where  K( aq,  aq,  aq,  aq)  =  r)(x\,  aq,  aq  =  Ar3,  aq,  aq).  We  now  adopt  the  convention  of  saying 
that  /(•)  is  the  integral  transform  of  the  function  /(•)  by  the  kernel  K(-).  The  factors  that 
determine  the  form  of  the  kernel  and  hence  the  type  of  transform  are  differential  operator  L 
and  the  coordinate  system,  which  in  turn  is  in  most  cases  determined  by  the  geometry  of  the 
boundary.  The  problems  we  are  interested  in  here,  is  in  mechanics  and  then  L  [■]  is,  in  general, 


1  We  assume  simple  boundary  conditions,  where  only  u  is  prescribed  at  the  boundary.  Other  types  of  boundary 
conditions,  involving  derivatives  of  u  can  be  treated  similarly. 
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the  Navier’s  operator: 


L[u]  —  (A  +  fi)A u  +  pV  x  V  x  u  —  pu  —  0  u  6  D,  (2  2) 

where  A,  p  and  p  are  material  constants  defined  later.  We  restrict  us  now  to  homogeneous 
materials  and  for  the  time  being  make  the  following  assumptions  on  boundary  conditions: 

1.  The  boundary  conditions  are  static,  so  dynamical  effects  can  be  ignored,  (u  =  0). 

2.  They  are  prescribed  in  Cartesian  coordinates,  then  L[-\  will  have  constant  coefficients, 
and  we  write  {x\,X2,x$)  =  (x,y,z),  . 

3.  Specifically,  let  the  boundary  conditions  be  u(x,y,  0)  =  f(x,y )  and  2  €  1R0+,  i.e.  the 
problem  of  a  semi-infinite  solid  or  a  halfspace  loaded  on  the  boundary. 

The  method  of  separation  of  variables  then  gives  us  that  Eq.  2.2  has  a  solution  u  =  exp(j(o;xa-  + 
u}yy)exp  —  which  is  valid  for  all  real  (ujX:  u>y).  From  which  we  obtain  the  kernel 

K  =  exp  (j(iOxx  +  (jJvy))  and  hence  recognize  the  integral  transform  as  being  the  inverse  Fourier 
transform  and  the  separation  constants  (uix,  ujy)  as  being  the  corresponding  spatial  frequencies. 

We  will  now  explore  some  the  properties  of  the  Fourier  transform  and  some  of  its  relatives. 


2.1  The  Fourier  Transform 

The  Fourier  transform  is  useful  when  when  the  DE  has  constant  coefficients  and  the  domain  of 
the  solution  extends  to  infinity  in  both  directions  (D  =  IR).  It  is  defined  as 

/+00 

f(X)e-’“*xdx.  (2.3) 

-OO 

Its  inverse  is, 

i  r+ 00 

f(x)  =  [/]  =  —/  f{ux)e^xdux.  (2.4) 

J  —  00 


Convolution  Theorem 


The  convolution  theorem  is  a  powerful  theorem  which  states  that  the  transform  of  a  convolution 
of  two  functions  is  the  product  of  the  transforms  of  each  function,  i.e. 


T 


/+00 

f(u)g(x  -  u)du 

-OO 


=  =  f(ux)g(u>x) 


(2.5) 
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Transfer  Functions:  This  theorem  is  very  useful  when  applied  in  conjunction  with  the  su¬ 
perposition  principle.  Convolution  can  be  viewed  as  a  expression  that  is  obtained  when  the 
principle  is  applied  to  a  shift  invariant  system  in  order  to  obtain  the  solution  for  a  general 
boundary  condition  from  a  solution  obtained  by  describing  the  boundary  conditions  as  an 
Dirac’s  impulse  function  6(x).  The  power  of  this  concept  lies  in  the  fact  that  one  only  has  to 
solve  the  problem  for  one  type  of  boundary  conditions,  the  impulse  function,  which  is  often 
relatively  simple.  A  general  solution  can  then  be  obtained  from  the  impulse  response  using- 
convolution,  which  becomes  a  simple  product  in  the  frequency  domain.  The  Fourier  transform 
of  the  impulse  response  is  defined  as  being  a  Transfer  Function ?  a  term  which  we  will  use 
frequently. 


Properties 


The  relation  between  the  behavior  of  a  function  and  its  Fourier  transform  is  in  many  ways 
interesting  and  several  theorems  exist  describing  that.  However,  here  we  will  only  detail  the 
relations  most  relevant  to  subsequent  analysis. 


Oddness  and  Evenness:  The  function  symmetry  properties  around  the  origin  have  some 
significance  as  for  real  functions  the  following  applies, 

If  f{x)  is  real  and  even  f(utx)  is  real  and  even 

If  fix)  is  real  and  odd  f(u>  x)  is  imaginary  and  odd 


Derivatives:  Differentiating  Eq  (2.4)  wrt.  x,  a  relation  for  the  Fourier  transform  of  derivatives 
is  obtained  as 


T 


dnf 
_  dxn  _ 


=  =  {jux)nJM 


(2.6) 


Scaling:  The  similarity  theorem  describes  how  scaling  of  original  function’s  variable  translates 
into  frequency  space,  and  it  states  that 

•H/(a*)l  =  "/(  —  )■  (2.7) 

a  a 

In  plain  english  this  means  that  a  function  that  is  ’flat’  has  a  Fourier  transform  that  is  ’sharp’ 
in  frequency  space,  and  reverse,  as  is  illustrated  in  Fig.  2-1. 

2For  spatial  systems  (not  dynamic)  the  impulse  response  is  frequently  called  point  spread  function  and  the 
Fourier  transform  of  that  the  modulation  transfer  function  [27]. 
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Spatial  domain  a=3 


Frequency  domain  a-3 


x 


Spatial  domain  a=//J 


jr 


w/(2n) 

Frequency  domain  a=l/3 


w/(2n) 


Figure  2-1:  Illustrating  the  effects  of  scaling  on  the  relationship  between  a  function  and  its  Fourier 
transform.  The  figures  on  the  left  are  the  scaled  functions  in  the  spatial  domain,  i.e.  f(ax)  where 
/(.'■')  =  exp(  —  irx2)  arid  the  figures  on  the  right  and  side  are  the  corresponding  Fourier  transforms,  which 
are  obtained  using,  3F(f)  =  exp(— 7ru;2).  The  Fourier  transform  is  then  scaled  in  the  frequency  domain 
according  to  Eq.  (2.6),  which  gives  us  that  F{f{ ax ))  =  |a|_1  f(u>x/  |a|)  =  |a|_1  exp(-7r(wz/a)2). 


2.2  The  2-D  Fourier  Transform 


The  2-D  Fourier  transform  is  a  generalization  of  the  1-D  version  above,  and  is  mainly  applicable 
to  partial  differential  equations  (PDE).  It  and  its  inverse  are  defined  as 


f{u;x,uy)  =  d  [/]  =  J I 

/<*,,)- *5!  l/l -£/£"/(■ 


+oo 


f(x,  y)e  ^u''xX+u,yy^,dxdy 

x,ujy)e^*x+^dujxdLOy. 


(2.8) 


Properties 


The  theorems  and  properties  delineated  for  the  1-D  Fourier  transform  above  can  be  generalized 
for  the  2-D  transform  and  are  in  brief: 

Oddness  and  Evenness: 
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Spatial  domain  a-3/2  b-2/3 


Frequency  domain  a=3/2  b=2/3 


Figure  2-2:  Illustrating  the  relationship  between  an  asymmetric  distribution  and  its  2D-Fourier  trans¬ 
form.  The  figure  on  the  left  is  the  distribution  in  the  spatial  domain,  i.e.  f(x,  y)  =  exp(— n(a2x2  +b2y2)) 
with  a  =  3/2  and  b  =  2/3.  The  figure  on  the  right  hand  side  is  the  corresponding  Fourier  transform, 
which  is  JFeP(/)  =  exp(-7r(oF/a2  +u >l/b2))/  |ai>|. 


If  f(.r,  y )  is  real  and  even 
If  f(x,y)  is  real  and  odd 


f(uix,ujy)  is  real  and  even 

is  imaginary  and  odd 


Derivatives: 


J~2D 


■  Qn  Qmy 

.  dxn  dym . 


=  (MnUuy)mr[f]  =  (M)n(M)m/(^.^) 


(2.9) 


Scaling:  The  two-dimensional  version  of  the  Similarity  theorem  states  that 

3~2D  [ f(ax,by )]  =  — ,  y-).  (2.10) 

|«o|  a  b 

Again,  this  means  that  a  function  which  is  concentrated  along  some  preferred  direction  has  a 
Fourier  transform,  which  preferred  direction  in  frequency  space,  is  generally  at  right  angles  to 
that  direction  (see  Fig.  2-2). 


2.3  The  Laplace  Transform 

Related  to  the  Fourier  transform  is  the  Laplace  transform,  and  it  is  applicable  when  the  domain 
D  of  the  solution  extends  to  infinity  in  only  one  direction,  i.e.  D  =  ffi°+  and  it  is  primarily 


3By  ‘odd’  we  mean  that  f(x,y)  has  the  same  sign  as  xy,  i.e.  positive  in  1st  and  3rd  quadrants  but  negative 
in  2nd  and  4th  quadrants. 
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used  on  dynamic  DE  when  the  effects  of  transients  are  included.  It  is  defined  as 


f  OO 

F(s)  =  £[/]=/  f(t)e~stdt 
Jo 

Its  inverse  is, 

1  rc+joo 

m  =  ^  lF}  =  ~r  F(s)estds 

JZ7T  J  c—joo 


(2.11) 


(2.12) 


Properties 


The  only  property  of  the  Laplace  transform  that  we  will  detail  here  is  how  derivatives  get 
transformed,  which  assuming  zero  initial  conditions  is: 


Derivatives: 


C 


d V 

dtn 


=  snIl[f}  =  snF(s). 


(2.13) 


2.4  Other  Transforms 

Finite  and  Discrete  Fourier  Transform 

One  of  the  properties  of  the  semi-infinite  halfspace  is  that  the  frequency  spectrum4  of  its 
solutions  is  continuous.  Conversely,  if  the  domain  D  is  finite,  say  x  €  [—L/2,  L/2],  the  spectrum 
becomes  discrete  as  then  the  solution  will  only  be  valid  for  discrete  values  of  u>x.  The  inversion 
formula  (Eq.  2.4)  therefore  becomes 

+  OO 

f(x)=  Y  f{n2n/L)eP^  (2.14) 

n=— oo 

and  the  limits  of  integration  in  Eq.  (2.3)  become  finite  ([—L/2,  L/2]). 

If  we,  additionally,  assume  that  we  have  only  available  sampled  values  ,  i.e.  discrete  in  x, 
the  transform  in  Eq.  (2.3)  becomes  a  finite  series  and  the  series  in  Eq.  (2.14)  also  becomes 
finite,  i.e. 

f(kds)  =  Y  f  (n/ds)  ej27r^  f(n/ds )  =  Y  /  (kds)  (2.15) 

n= 0  0 

where  ds  is  the  sampling  period. 

One  can  also  derive  derivative,  scaling,  convolution  theorems,  similar  to  Eqs.  (2.6),  (2.7), 
and  (2.5)  for  the  finite  and  discrete  Fourier  transforms  and  also  a  discrete  and  finite  2D- 

4Spectrum  is  defined  as  being  the  set  of  frequencies  for  which  a  solution  is  valid. 
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Fourier  transform  [5].  The  derivations  are  straightforward  and  the  behavior  is  similar,  as  in  the 
continuous  case,  except  for  the  convolution  theorem,  where  the  effects  of  aliasing  appear.  Care 
has  to  be  taken  and  one  has  to  be  aware  of  the  effects  of  aliasing  when  results  given  here  are 
being  implemented  or  simulated  with  a  computer. 


Cylindrical  and  Spherical  Geometries 


When  the  geometry  is  such  that  the  equilibrium  equations  are  most  suitably  expressed  in 
cylindrical  or  spherical  coordinates,  the  same  solution  procedure  can  be  applied. 

Cylindrical:  In  the  case  of  cylindrical  geometry,  the  solution  can  be  expressed  as  a  combination 
of  a  Fourier  and  a  finite  Fourier  transform  and  a  series  involving  Bessel  functions,  e.g. 

00  /•+ 00 

u(r,e,z)  =  J2  /  T(r,n,LUz)Jn(LU,r)f(n,u}z)ejn6+J0JzZduz  (2.16) 

n= 0J-°° 

where  Jn(-)  is  a  Bessel  function  of  the  first  kind,  of  order  n.  /(•)  is  determined  by  the  boundary 
conditions,  and  T(-)  can  be  viewed  as  having  the  role  of  a  transfer  function. 

Spherical:  When  the  equilibrium  equations  are  expressed  in  spherical  coordinates  (r,  0 ,  <p),  the 
solution  can  be  expressed  as  a  series  involving  Legendre  polynomials,  e.g. 

OO  OC 

«(r.M)  =  £  £  nr<n,m)Pn{cosOym*Kn,m)  (2.17) 

n=0  m=  —  oc 


where 


pm 

*71 


2  \m 


W  =  2^!(1  -  "  > 


/2dm+n(ti2-iy 
(lnm+n 


is  the  Rodriguez’s  formula  for  the  associated  Legendre  polynomials  and,  as  before,  /(n,m),  is 
determined  by  the  boundary  conditions  and  T(-)  can  similarly  be  viewed  as  a  transfer  function. 


3.  Static  Indentation  Problems 

In  this  section  we  will  demonstrate  horv  the  mathematical  tools  developed  in  the  last  section  can 
be  applied  to  obtain  solutions  for  static  contact  problems.  We  start  by  giving  examples  of  simple 
system,  but  gradually  make  the  models  more  detailed.  By  this,  we  hope  to  give  the  reader  a 
feel  for  the  role  of  each  term  in  the  solution  and  that  in  turn  will  make  the  final  solution  less 
obscure  and  more  believable.  Having  achieved  that,  we  believe  that  we  have  demonstrated  the 
power  of  transfer  functions  and  integral  transforms  to  provide  understanding  of  the  qualitative 
behavior  of  the  solutions  to  contact  problems  by  viewing  them  as  the  output  from  a  system 
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Figure  2-3:  Shows  a  schematic  diagram  of  the  Winkler  elastic  foundation.  The  foundation  shown,  is 
being  indented  with  a  rigid  object  with  a  total  load  force  Fz,  resulting  in  a  load  distribution  fz(x).  The 
material  is  modeled  as  a  collection  of  non-interconnected  springs.  The  loading  at  each  point  on  the 
surface  contributes  therefore  only  to  the  displacement  at  that  location  and  does  not  extend  beyond  the 
contact  region. 


which  has  the  boundary  conditions  as  inputs. 


3.1  Winkler  Foundation 

A  very  simple  way  of  modeling  a  compliant  material  is  the  Winkler  foundation,  where  the 
material  is  modeled  as  a  collection  of  non-interconnected  identical  springs,  (see  Fig.  2-3)  and 
the  effects  of  shear  are  ignored.  The  force  displacement  relations  we  obtain  are  simply 

ID  :  uz(x)  =  2D  :  uz(x,y)  =  ^-fz(x,y)  (2.18) 

'yjU) 

where  uz  is  the  displacement  in  the  2  direction,  kw  the  stiffness  per  unit  length  (ID)  or  area 
(2D)  of  the  foundation,  and  fz  is  the  normal  loading  per  unit  length  (ID)  or  area  (2D).  The 
relations  in  transfer  function  form  are  identical. 
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3.2  Strings  and  Beams  under  Tension  on  a  Winkler  Foundation 


The  Winkler  foundation  model  has  the  advantage  of  being  simple,  but  the  fact  that  any  material 
deforms  continuously,  and  that  the  deformations  are  therefore  not  independent  is  neglected.  The 
springs  should  therefore  somehow  be  interconnected.  One  way  to  include  this  effect  is  to  model 
the  material  as  a  string  under  tension,  S,  lying  on  a  Winkler  foundation.  We  assume  that  the 
nature  of  the  string  is  such  that  it  can  only  resist  load  along  its  length5.  The  energy  stored 
is  related  to  the  elongation  of  the  string  which  is  related  to  the  first  derivative  of  the  surface 
deformation  and  the  principle  of  virtual  displacements  gives  us  the  equilibrium  equation  [61]: 

-  +  kwuz  =  fz  (2.19) 

which  in  transfer  function  form  is 


uz{lox) 


Suj2  +  kx 


-fz(ux), 


(2.20) 


where  S  is  the  string  tension.  The  above  can  be  generalized  further  to  include  higher  order 
derivatives  of  uz  and  if  we  model  the  medium  as  a  beam  under  tension  S  with  bending  stiffness 
El,  resting  on  a  Winkler  foundation;  then  using  similar  methods  we  obtain  the  transfer  function 
relation 


u  z  (T'.f ) 


1 

Elu^i  +  Slu2  +  kw 


fz(vx)- 


(2.21) 


3.3  Membranes  and  Plates  under  Tension  on  a  Winkler  Foundation 


The  above  models  can  readily  be  generalized  to  describe  the  deformations  of  a  surface  sub¬ 
jected  to  two-dimensional  loading  where  a  string  corresponds  to  a  stretched  membrane,  and 
a  beam  corresponds  to  a  plate  under  tension  S ,  having  flexural  stiffness  N.  In  that  case  the 
corresponding  equilibrium  equations  are  PDEs  in  x  and  y  and  hence  we  use  Eq.  (2.9)  to  obtain 
the  solution  in  a  transfer  function  form. 


1 


Uy)  —  _|_  £  fz(u,xil^y)i 


and 


(Wl,  ix)y  ) 


ArQ4  +  SO,2  +  fc, 


'  f )(wi!  w y)i 


where  Q  =  ^ 


(2.22) 

(2.23) 


“The  string  is  therefore  equivalent  to  having  the  load  points  connected  to  its  nearest  neighbor  with  a  spring 
in  tension. 
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Remark:  That  we  express  the  above  results  in  terms  of  transfer  function  enables  us  to  interpret 
the  relation  between  the  normal  loading  and  displacement  as  spatial  low-pass  filtering.  The 
bandwidth  of  the  filter  is  then  prescribed  by  the  stiffness  and  the  tension;  higher  stiffness  or 
tension  translating  to  lower  bandwidth  and  therefore  smoother  deformations. 

3.4  Infinite  Half-space 

In  the  above  we  have  only  obtained  expressions  for  the  surface  deformations  of  the  medium, 
subjected  to  normal  load.  Since,  we  will  in  the  chapters  to  follow  need  a  model  that  can  be  used 
for  the  analysis  of  the  response  of  sub-surface  sensors  we  need  to  extend  our  model  to  include 
sub-surface  deformations  and  also  take  into  account  tangential  or  shear  loads.  We  therefore 
model  the  medium  now  as  being  homogeneous  and  extending  to  infinity6  in  the  2  direction. 

In  the  case  of  a  continuum  model,  the  solution  are  not  obtained  as  directly  and  easily  as  for 
the  simple  foundation  problem  treated  above,  and  we  resort  to  the  following  procedure7: 


Solution  Procedure 

1.  The  starting  point  is  the  equilibrium  equations  and  the  stress-strain  relations  (generalized 
Hooke’s  law). 

2.  We  then  combine  the  two  sets  of  equations  to  eliminate  the  stresses  and  obtain  cou¬ 
pled  second  order  PDE’s  in  terms  of  the  displacements  only,  effectively  writing  out  the 
transformed  Navier’s  equation. 

3.  We  then  transform  the  PDEs  using  the  Fourier  transform  with  respect  to  the  appropriate 
variables  (in  general  x,y).  The  PDEs  are  then8  transformed  into  a  system  of  coupled 
ODEs  in  2. 

4.  The  ODE  system  is  solved  using  standard  methods  with  prescribed  surface  loading  or 
displacements  as  boundary  conditions. 


6 The  solution  method  used  here  can  also  be  extended  to  aelotropic  (layered)  and  finite  thickness  materials, 
but  the  solutions  quickly  become  quite  complicated  and  opaque  and  further  discussion  will  therefore  not  be  given 
here. 

'  For  the  sake  of  simplicity,  we  assume  here  that  the  solution  domain  is  the  semi-infinite  half-space,  and 
therefore  most  easily  described  using  Cartesian  coordinates.  Solution  methods  for  other  domains  can  be  obtained 
similarly,  e.g,  as  shown  in  section  (2.4). 

8Steps  2  and  3  can  be  interchanged. 
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Plane  Strain  (Line  Load) 


Example:  We  will  show  how  the  solution  for  an  elastic  half-space  loaded  over  a  narrow  strip,  so  that 
the  loading  can  be  considered  a  distribution  of  line  loads.  We  choose  our  coordinate  system  so  that  the 
boundary  surface  coincides  with  the  xy  plane  and  the  2-axis  is  directed  into  the  solid.  The  tractions 
(normal  and  tangential)  are  assumed  not  to  vary  along  the  y-dimension  and  hence  be  functions  of  '$ 
only.  We  shall  assume  that  that  under  these  conditions,  a  state  of  plane  strain  ( eyy  —  0)  is  produced  in 
the  half-space. 

Stepping  through  the  procedure  introduced  above  we  obtain: 


Step  1:  Under  the  simplifying  plane  strain  assumption, 


do 


dx 

doza 

dx 


xx  ^  dox 


+ 


dz 

doz 


dz 


the  equilibrium  equations  become: 
=  0 

=  0 


and  the  stress-strain  relationships: 


(2.24) 


Cn  —  (A  +  2  /.i)sXx  +  As--,  crzz  —  (A  +  2y)ezz  +  AsIX,  oxz  —  2/jexz  (2.25) 


where 


A  = 


E  v 

(l  +  i/)(l-2i/) 


and  fi  = 


E 

2(l  +  i/) 


(2.26) 


are  the  Lame  constants  (bulk  and  shear  modulus)  with  E  and  v  denoting  Young’s  modulus  and  the 
Poisson  ratio,  respectively. 


_  dux 

&xx  —  7  ^zz 

OX 


duz 
dz  ’ 


£xz 


1  (  dux  duz  \ 

2  {~dT  +  ~fa) 


(2.27) 


Step  2+3  Carrying  steps  2  and  3  out  simultaneously  we  get  the  system  of  equations 

j(p'2  -  l)Dujx 
j(/32-l)Dwx  02D2-ul 
d 


£ 

O 

£ 

H 

o 

(2.28) 


where  (52  =  (A  +  2//)//x  and  D  = 


dz' 


Step  4:  The  determinant  of  the  system  in  Eq.  (2.28)  is  A  =  02(D2  —  u>l)2,  and  a  non-trivial  solution 
will  therefore  be  of  the  form 

ux  =  (Ax  +  Bxz)e^^z 

(2.29) 

uz  —  (.4;  +  B2z) 


which  has  4  unknown  coefficients.  Two  of  them  can  be  eliminated  by  substituting  Eq.  (2.29)  back  into 
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Eq.  (2.28)  which  gives: 


Bx  =  -jsgn(ux)Bz,  A-  =  jsgn(ujx)A: 


+  p2  +  1B 
+  pTTiB> 


where  sgn  is  the  sign  function  with  sgn(O)  —  0.  The  two  remaining  coefficients  are  decided  by  the 
boundary  conditions,  which  we  assume  are  in  the  form  of  prescribed  tractions  on  the  boundary  surface 
(z  =0).  The  equations  relating  the  coefficients  to  the  tractions  are  obtained  by  substituting  Eq.  (2.29) 
into  Eqs.  (2.25)  and  (2.27)  and  identifying  azz(z  =  0)  as  the  normal  traction  and  ax,(z  =  0)  as  the 
tangential  traction.  We  then  get  the  relations 

2/x ((1  -  2t/)sgn (ux)Bz  -  |wx|  Ax)  =  fx 

2/J  ((2  -  2v)Bz  -  juxAx)  =  /,. 


Solving  for  Ax  and  Bz  and  substituting  back  into  Eq.  (2.29)  gives 


«,  =  ■  P-2pklUe-|„|7  +  1(1  -2--KI--)e-t,,|,/  =  T  j  +  T  f 

2//  \oJx\  2ft  u >x 


Uz  =  --j 


+  A.  i2.jzlv  +  J.^1  Zle-WA--  f  =T  f  f  +T  f  f 

C  JlT  ||  C  Jz  —  1  uz  fx  J  X  ^  -L  U~fz  J  Zi 

.{X  L Ox  I  1 


(2.30) 


where  TUjjj:  and  TUtfj  are  the  transfer  function  for  the  displacements  ux  and  uz,  resulting  from  the 
tangential  load  fx,  and  similarly  TUrft  and  TUrfr  are  the  transfer  functions  for  the  same  displacements 
displacements,  resulting  from  the  normal  load  fz. 


We  express  this  solution  compactly  using  a  Transfer  Function  Matrix  (TFM),  i.e. 


u  =  Tuff ,  (2.31) 

where  u  =  [ux,  u2]r,  f  =  [fx,  fz]T ,  and 


"  2  -  2u  —  \lox\z 

.1  —  2 V  —  |tux|  2  " 

Tux  fx 

Tux  fz 

e 

ksl 

^  W.T 

1 - 

b" 

Tuz  fz  _ 

2  fi 

1  —  2  v  +  \lux  1  2 

J 

2  —  2v  +  \ux\  2 

W. X 

\^x ! 

Assuming  that  fz(x)  is  an  even  function,  we  observe  that  the  resulting  ux(lox)  is  imaginary 
and  odd,  indicating  that  ux(x)  is  odd  and  similarly  we  conclude  that  the  u~(x)  resulting  from 
an  even  fz  is  also  even. 

Symmetry:  The  terms  in  the  TFMs  appear  similar  but  there  are  subtle  differences  so  they 
are  not  completely  symmetric  in  ux  and  2.  This  can  be  expected,  mainly  because  u jx  is  a 


28 


frequency  variable  but  z  is  a  spatial  variable.  Also  since  the  domains  of  each  are  different 
(ujx  G  [-00,  +00],  but  z  G  [0,  +00]). 


The  strains  and  stresses  can  then  be  obtained  using  Eqs.  (2.25)  and  (2.27)  respectively, 
giving: 


and 


&XX 

j sgn(cjI)(2  -  2v  +  \ux\  z) 

1 - 

1 

CM 

i 

T— 1 

jsgn(u>x)(2iy  -  \u>x\z) 

1  —  2is  +  \ujx\  z 

SZz 

2p 

&XZ 

1  -  |wx|  z 

~jiOxZ 

fx 

fz 


&XX 

-jsgn(u;x)(2  -  |wx|  z) 

1  -  \ux\ Z 

Ozz 

=  e-l^t~ 

~juxz 

1  +  \wx  \  z 

&XZ 

1  l^'xl  Z 

j  vxz 

fx 

fz 


(2.32) 


(2.33) 


Mean  normal  Stress:  A  stress  combination  that  becomes  particularly  important  for  incom¬ 
pressible  materials  is  the  mean  normal  stress,  defined  as 

A  1  ,  ,  x 

Pn  —  g  ®yy  '  ®  zz  J  • 


Its  importance  lies  in  that,  for  incompressible  materials,  it  does  not  contribute  to  any  strains, 
(an  incompressible  body  subjected  to  uniform  pressure  does  not  deform  at  all).  This  means  that 
strains  cannot  be  used  solely  as  a  sensor  signal,  if  one  wants  to  use  incompressible  materials 
as  a  cover  for  tactile  sensors,  as  then  the  sensors  would  not  be  able  to  sense  the  mean  normal 
stress  component  of  the  stress  signal. 

For  the  line-load  case  we  get  that  pn  =  (crxx  4-  <jzz)  since  then  ayy  =  v(crxx  +  azz)  and 

using  Eq.  (2.33)  we  get  the  TFM  for  pn  as 


Tp„f  — 


2(1  +  v)e~^z 


-j  sgn(u4), 


1] 


Other  Structures:  Comparing  the  solution  for  the  z-displacements  (u2)  on  the  surface  (z  = 
0)  of  the  halfspace  subjected  to  distributed  normal  load  to  those  obtained  for  the  Winkler 
foundation  and  the  string  we  observe  that: 


M  Winkler  ^  fz 


^Halfspace  ^ 


^String  OC  of' 


vt 
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Filtering  characteristics  of  transfer  functions 


Figure  2-4:  Illustrating  the  shape  of  typical  terms  in  the  transfer  functions.  The  dash-dot  curve  (i) 
plots  T  =  exp(—  |oij;|  z)  =  TPnft  which  has  lowpass  characteristics.  The  dashed  line  (ii)  plots  T  oc 
|cax|  c exp(—  \u)x\ z)  =  Trtifz  which  has  bandpass  characteristics.  The  solid  line  (iii)  plots  these  two 
added  or  T  —  (1  +  |u;.r|  c)exp(—  \uix\  z)  =  which  has  increased  bandwidth  over  (i) 

The  relative  degree  of  the  transfer  function,  tells  us  how  localized  the  resulting  deformations 
will  be.  The  Winkler  foundation  deforms  only  directly  under  the  load  and  that  is  reflected 
in  a  flat  spectrum,  indicating  a  very  sharp  impulse  response.  The  surface  displacements  of 
the  halfspace  spread  away  from  the  loading  area,  because  the  transfer  function  goes  to  zero 
for  higher  frequencies  the  medium  must  therefore  deform  gradually.  The  string  has  an  even 
sharper  transfer  function  so  its  displacements  must  spread  out  even  further.  The  semi-infinite 
halfspace  therefore  behaves,  in  terms  of  its  surface  displacement  profile  response  to  distributed 
normal  loading,  as  something  that  is  between  a  Winkler  foundation  and  a  string  under  tension. 

Spatial  Filtering:  The  exponential  term  represents  spatial  lowpass  filtering  which  bandwidth 
decreases  exponentially  with  depth.  For  the  linear  strains  and  stresses,  two  main  variants  of 
terms  can  be  identified  in  Eqs  (2.32)  and  (2.33):  i)  e~^z,  ( pn )  ii)  ^\u/x\e~^z,  (ezz).  The 
shapes  of  each  are  plotted  in  Fig.  2-4.  There  we  see  that  (i)  corresponds  to  a  low  pass  filter 
but,  (ii)  a  bandpass  filter.  The  combination  of  (i)  and  (ii)  (1  +  2  |tnx|)e~lWxl~  is  also  plotted  as 
a  solid  curve  and  it  is  seen  that  it  covers  both  the  low  and  high  frequency  bands. 
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pn(x)  [Normalized] 


(i)  T  =  expC-lcoJz) 


A 


(ii)  T  =  dco(.lexp(  — lcovlz) 


x 


X 


Figure  2-5:  Shows  the  response  in  the  spatial  domain,  of  each  of  the  transfer  functions,  (i)-(iii)  to  uniform 
normal  loading  in  x  G  [—1/2, 1/2].  In  the  leftmost  figure  is  the  mean  normal  pressure pn  —  —  (<rII+crZJ)/2 
is  plotted.  It  has  the  transfer  function  T  =  exp(—  |cax|  z).  In  the  middle  figure  the  ezz  normal  strain  is 
plotted,  it  has  transfer  function  T  =  z|u>x|exp(—  |ujx|  z).  In  the  rightmost  figure  the  cjz-  normal  stress 
is  plotted.  It  has  transfer  function  T  —  (1  +  z  |u>x|)exp(—  |o>x|  z). 

For  the  spatial  domain,  the  effects  of  the  different  filtering  functions  are  shown  in  Fig.  2-5 
and  there  we  see  that  (i)  and  (iii)  are  similar,  but  (i)  is  more  flat.  The  (ii)  terms  on  the  other 
hand  show  typical  bandpass  characteristics  such  as  sign  changes  near  the  edge  of  the  contact 
region  (x  =  ±1/2). 


Sampling  period  vs.  Sensor  Depth:  We  also  note  the  inverse  relationship  between  z  and  ux 
which  becomes  quite  important  when  one  considers  the  consequences  of  using  discrete  sensors 
and  hence  the  effects  of  sampling.  It  turns  out  that  it  has  an  equivalent  effect  on  the  bandwidth 
of  the  transfer  function,  to  either  place  the  sensors  further  apart  or  further  beneath  the  surface. 
This  indicates  a  tradeoff  between  sensor  depth  and  sensor  spacing  which  becomes  important 
when  sensor  arrays  are  designed. 


General  3D  Solution 

Exactly  the  same  procedure  can  be  applied  when  the  plane  strain  assumption  is  dropped  and 
general  loading  conditions  are  considered.  That  the  tractions  vary  also  along  the  ^-dimension 
is  accommodated  by  using  the  2D  Fourier  transform  defined  in  Eq.  (2.8).  The  solution  can 
similarly  be  presented  using  TFMs  with  full  dimensions  of  the  displacement,  traction  and 
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stress/strain  vectors,  i.e. 


£  —  [£xxi  £yy>  &xy>  £xzi  £yx J 

cr  [^ii,  ®yyi  0”-2,  @xyi  r^yx]  • 

The  details  of  the  TFMs  will  be  omitted  here  but  the  full  solution  can  be  found  in  chapter  3. 
All  the  comments  made  above  regarding  the  plane  strain  solutions  apply  directly  to  the  general 
3D  solution.  Although  one  has  to  keep  in  mind,  when  transferring  the  concept  of  bandwidth 
over  to  the  2D  transfer  functions,  that  the  frequency  space  is  two-dimensional.  The  bandwidth 
of  transfer  functions  will  therefore  in  general  be  dependent  on  the  direction  also. 


U  [U'x i  ^yi  ^z\  i 
f  =  [/*,  fy,  fz]T  , 


4.  Dynamic  Problems 


The  problems  we  have  discussed  so  far  have  assumed  static  loading,  i.e.  we  have  ignored  both 
the  effects  of  rate  dependent  material  laws  (viscosity)  and  also  the  effects  of  inertia,  but  we 
can  treat  these  two  phenomena  separately.  The  effects  of  inertia  and  wave  mechanics  will 
be  discussed  separately  for  each  structure  and  we  will  also  assess  when  inertial  effects  can  be 
ignored.  That  will  simplify  our  solution  significantly,  especially  in  the  case  of  a  semi-infinite 
solid.  The  viscoelastic  behavior  is  modeled  by  modifying  the  material  laws  and  making  them 
rate  dependent. 

The  effects  of  inertia  p  are  included  by  adding  d’Alemberts  force  f  \  =  puz  to  the  right 
hand  side  of  the  equilibrium  equations.  Using  Laplace  transform,  the  double  time  derivative 
can  be  replaced  with  a  multiplication  by  s 2  and  the  only  addition  to  the  static  solution  for  the 
foundation  models  is  a  extra  term  of  ps 2  in  the  denominator  of  the  transfer  function.  We  will 
discuss  the  effects  of  inertia  on  the  solution  to  the  structures  analyzed  in  sections  3. 1-3.3  in 
the  sections  4. 1-4.3.  The  inertia  terms  have  more  complex  effects  on  the  semi-infinite  halfspace 
solution  and  we  will  take  a  look  at  those  in  section  4.4.  How  viscoelasticity  can  be  included  in 
the  transfer  function  formulation  will  finally  be  discussed  in  section  4.5. 


4.1  Winkler  Foundation 


For  this  simple  case  the  solution  now  becomes: 


uz 


1  1 

p  <jj\ ]  +  s2 


fz. 


(2.34) 
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where  p  is  the  density,  and  ojq  —  kw/p  is  the  resonance  frequency  of  the  foundation. 


4.2  The  String  and  the  Membrane 


For  the  string  the  solution  is: 


and  for  the  membrane: 


uz 


1  1 

p  c%tol  +  s2 


h 


Uz 


1  1 

p  +  s2 


h 


where  ci  —  \/T j p  is  the  wrave  velocity  as  can  be  seen  if  we  assume  the  loading  /.  is  a  point 
load  at  x  =  0  but  then  the  Laplace  transform  of  the  solution  for  the  string  becomes 


uz(x,s)  =  - 


C2  ■ -^h(s) 

— e  c2  - 

2s  T 


We  observe  that  viewed  from  the  point  of  loading  the  response  of  the  string  is  that  of  a  pure 
damper.  This  is  a  consequence  of  the  string  having  no  transverse  or  bending  stiffness,  and 
hence  all  the  power  injected  into  the  system  is  transferred  away  from  the  source  (loading  point) 
and  appears  lost.  The  exponential  term  indicates  that  the  response  away  from  the  origin  is  a 
time-delayed  copy  of  that  at  the  loading  point,  but  delayed  by  =  |.?:|  jc2  which  is  according 
to  to  the  well  known  result  of  d’Alembert  [61],  namely  that  the  solution  is  a  wave  traveling  in 
both  directions  with  velocity  C2  away  from  the  origin. 


For  the  membrane,  we  obtain  for  a  time-varying  point-load  at  the  origin  that 


_  ,  .  27r/,n(s)  r_  ,rs. 

uz(r,  s)  =  —  J^Ao  - 
T  c2 


where  A'o(-)  is  a  modified  Bessel  function  of  the  second  kind.  Using  asymptotic  expansions  for 
the  Bessel  function  we  obtain  that  the  steady-state  response  to  a  sinusoidal  excitation  at  the 
origin  for  points  that  are  far  away  from  the  origin  is 


u.(r,t)  m  j(c2  w/4)  sin(u;f), 

T  V  irruj 


where  Fzo  is  the  total  load  and  u  is  the  excitation  frequency.  This  corresponds  to  a  traveling- 
wave,  with  a  phase  delay  of  ir/4.  and  with  amplitude  that  decreases  as  A-1/2  where  A  is  the 
distance  from  the  origin  measured  in  wavelengths  A2  =  ttco/’X. 
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4.3  String  or  Membrane  on  a  Winkler  Foundation 


If  the  string  or  the  membrane  is  supported  by  a  Winkler  foundation  the  transfer  function 
solution  becomes 


uz 


1  1 
P  + 52 


h 


for  the  string  and 

-  1  1  1 

p  +  Wq  +  s2  z 

for  the  membrane.  For  input  frequencies  less  than  the  resonance  frequency  of  the  Winkler 
foundation,  (|s|  <  uq)  the  response  is  that  of  a  pure  spring  but  for  higher  frequencies  it  behaves 
in  part  like  a  string  and  some  of  the  power  of  those  frequencies  gets  radiated  away. 


Neglecting  Inertial  Effects:  It  seems  reasonable  that  for  slow  motion,  the  effects  of  inertia 
could  be  neglected  as  they  would  be  overwhelmed  by  the  static  and  viscous  effects.  By  inspecting 
the  solution  we  see  that  the  effects  of  inertia  can  be  ignored  if: 

,  c2 

ojq  3>  u!  and  —  3>  u> 

L 

where  u.’  is  the  (maximum)  input  frequency  and  L  is  a  (maximum)  characteristic  length,  i.e. 
scale,  of  the  structure  of  interest.  The  first  condition  indicates  that  the  input  frequency  should 
not  excite  the  resonance  frequency  of  the  Winkler  foundation  and  the  second  condition  indicates 
that  the  wavelength  of  the  wave  excited  by  the  highest  input  frequency  of  the  input  should  be 
large  compared  to  the  scale  of  the  structure. 


Beams  and  Plates:  Solutions  for  dynamic  loading  of  beams  and  plates  can  be  obtained  using 
the  same  method  as  above.  No  analysis  of  their  behavior  will  be  given  here. 


4.4  Infinite  Half-space 

In  this  section  we  will  study  the  dynamic  behavior  of  the  infinite  halfspace.  First  we  will 
analyze  the  mechanics  of  the  different  types  of  plane  waves,  then  how  the  same  types  of  waves 
appear  from  a  point  source.  We  will  finally  analyze  the  dynamic  response  of  the  elastic  halfspace 
subjected  to  a  dynamic  line-load,  and  based  on  that  give  conditions  under  which  we  can  neglect 
the  effects  of  inertia. 
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Wave  Mechanics  of  Plane  Waves 


At  a  sufficient  distance  from  their  source,  all  waves  can  be  considered  plane  waves  and  it  is 
informative  to  consider  the  propagation  of  plane  waves.  A  plane  harmonic  displacement  wave 
can  be  represented  as 

u  =  Adexp  (jk(r  ■  p  —  ct))  (2.35) 


where  u  is  the  displacement  vector,  d  is  a  unit  vector  defining  the  direction  of  motion  of  the 
particle,  p  is  a  unit  vector  defining  the  direction  of  propagation,  c  is  the  wave  velocity,  k  is  the 
wave  number  and  A  is  independent  of  the  position  vector  r  and  time  t.  When  Eq.  (2.35)  is 
substituted  into  Navier’s  Eq.  (2.2)  we  find  that  two  types  of  plane  waves  are  possible,  when 
either 


(i)  d  =  ±p  and  c  =  c„  = 


A  A  2  p. 


1/2 


or  (ii)  d  ■  p  =  0  and  c  =  csh  = 


l± 

P  J 


1  1/2 


(2.36) 


In  case  (i)  the  motion  of  the  particles  is  parallel  to  the  direction  of  propagation  and  the  wave  is 
a  longitudinal  wave.  It  can  be  shown  to  involve  only  volume  dilation  and  no  shear  or  rotation. 
In  case  (ii).  on  the  other  hand,  the  motion  of  the  particles  is  perpendicular  to  the  direction  of 
propagation,  involving  only  shear,  hence  called  a  shear  wave. 


A  third  type  of  waves  occurs  close  to  the  surface  of  the  solid  as  a  consequence  of  an  interac¬ 
tion  of  both  the  dilational  and  the  shear  wave  with  the  free  surface.  They  are  called  Rayleigh 
waves  and  a  solution  for  them  can  be  found  by  substituting 

u  =  AdPek^rP-CRt)-apZ^  +  dshek^{r'P-CRt)-a^  (2.37) 


into  Eq.  (2.2)  and  also  use  the  condition  that  the  surface  is  stress-free,  i.e.  ozz  =  aTZ  =  0  when 
c  =  0.  It  can  be  shown  that  these  conditions  lead  to  an  equation  for  the  ratio  c\R/ash 

({cR/csh)2  -  1/2)  -  {cR/ csh)2  \J \cR/ csh)2  -  P~'2\J{cR/csh)2  -  1  (2.38) 

that  can  be  shown  to  give  exactly  one  solution  that  can  be  interpreted  as  a  wave  [81].  It  turns 
out  that  —  €  [0.8740,  0.9553]  for  v  £  [0,  0.5]  and  we  further  get  that  ap  =  [1  —  c^/cp]1/2, 

csh 

ash  —  [l  —  c2r/c2si j]1  •  The  nature  of  Rayleigh  waves  is  that  they  propagate  essentially  only 
on  the  surface,  as  they  attenuate  exponentially  with  depth.  Their  geometric  attenuation  is 
therefore  lower  than  that  of  the  dilational  and  shear  waves  and  the}’'  will  die  out  slower  with 
distance  from  the  load  source. 


Geometric  attenuation:  When  a  sinusoidal  point  load  is  applied  at  the  origin  P  =  Pq  cos (ut) 
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the  radial  and  transverse  displacements  can  be  written  [36] 


UR  =  ^zj^QR(Q,P)cos(Lo't  -  kpR) 

p 

ue  =  o  °pQe(^,/?)sin(u;t  -  kshR) 
2-krR 


where  kp  =  u/cp  and  ksh  =  u>/cSh,  but  Qr(-)  and  ©#(•)  depend  on  the  directional  angle 
0  =  cos ~1(z/R)  and  (3  =  ^Jcp/csh-  These  equations  are  valid  everywhere,  except  at  the  surface 
0  =  7t/2.  We  observe  that  the  radial  displacement  ur  is  due  to  the  dilational  wave  only, 
and  similarly  the  tangential  displacement  is  due  to  the  shear  wave  only.  At  the  surface  the 
displacements  are  due  to  the  Rayleigh  wave  and  are  given  by 
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where  kR  =  uj/cr  and  Fr(-)  and  Fz(-)  depend  on  Poisson’s  ratio  [51].  Similar  expressions  can 
be  obtained  in  the  case  of  a  line-load.  The  asymptotic  behavior  of  this  solution  is  determined 
by  the  geometric  attenuation,  which  is  the  result  of  the  increasing  area  of  the  wavefront  with 
distance  from  the  source.  The  geometric  attenuation  depends  on  how  many  dimensions  the 
wave  propagates  in,  and  we  therefore  have  different  attenuation  for  Rayleigh  waves  than  for 
dilational  and  shear  wraves.  The  attenuation  also  depends  on  the  source  being  a  line-  or  a 
point-load.  We  have  that  the  amplitude  .4  behaves  with  distance  from  the  source  A  (measured 
in  wavelengths)  as  follows  [36,  51]: 


Line- load 

Dilational  and  shear  waves:  .4  oc  — r-^r 
Rayleigh  waves:  .4  oc  1. 


Point-load 

Dilational  and  shear  waves:  A  oc  — 

A 

Rayleigh  waves:  A  oc  y~[y^r- 

Hence,  Rayleigh  waves  in  the  halfspace  behave  like  shear  waves  in  a  string  for  a  line-load 
source,  and  as  shear  waves  in  a  membrane  for  a  point-load  source. 
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Transfer  Function  Solution 


The  dynamic  response  of  the  semi-infinite  halfspace  is  considerably  more  complicated  than  that 
of  the  structures  in  the  previous  examples,  mainly  because  displacements  are  now  allowed  to 
vary  also  in  the  ^-dimension.  The  solution  procedure  is  however  the  same  as  used  for  the  static 
problems,  detailed  in  section  3.4,  and  as  above,  the  effects  of  inertia  are  included  by  adding 
d’Alemberts  force  /a{  =  pUi  i  =  x,  y,  z  to  the  right  hand  side  of  the  equilibrium  equations.  For 
the  example  in  section  (3.4)  this  leads  to  the  modified  system  of  equations 
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which  has  the  determinant  A  =  (32(Dz  —  n\){D\  -  n 2),  where  ni  =  (cu2  +  a^2m-s2)1/'2  anfl 
no  =  ( u. >1  +  ^s2)1/2.  A  non-trivial  solution  will  therefore  be  of  the  form 


um  =  Axe~niZ  +  Bxe~n2Z 
uz  =  Aze~niZ  +  Bze~n2Z 


(2.40) 


A  complete  solution  of  both  ux  and  uz  displacements  resulting  from  a  combined  tangential,  fx, 
and  normal  load,  fz  can  be  obtained  using  the  same  method  as  in  section  3.4  as  well  as  the 
general  3D-solution  including  uy  and  fy.  These  can  then  be  used  to  derive  the  stress  and  strain 
relations  as  needed.  We  will  leave  out  the  details  and  only  look  at  the  normal  displacement 
response  for  a  distributed  line-load  and  also  for  a  general  distributed  load,  since  an  analysis 
of  their  dynamical  behavior  will  directly  carry  over  to  the  other  displacements,  stresses  and 
strains. 


The  solution  for  the  normal  displacement  uz  resulting  from  a  line-load  is  in  transfer  function 
form, 


uz (ujx .  s)  — 


n  1 


■niz  _  ,2 p-n2z 


2/^  (WI  +  2^2)2  -  wlnln2] 


f z{^xi  S). 


(2.41) 


The  point-load  solution  is  almost  identical,  one  only  has  to  replace  uj2  with  Q2  —  uj2  +  cu2  in 
the  transfer  function  and  change  the  arguments  in  uz  and  /,  to  (cjx,cjy,s). 


Comparing  Eqs.  (2.35)  and  (2.36)  to  the  transfer  function  in  Eq.  (2.41)  we  see  that  the 
exponents  n*  and  n-i  correspond  to  the  dilational  and  shear  waves  respectively.  The  Rayleigh 
wave  on  the  other  hand,  corresponds  to  the  pole  of  the  denominator  [51].  The  exponential  terms 
e~niZ  and  e~n2Z  correspond  to  the  interaction  of  a  traveling  wave  with  the  effects  of  increasing 
low-pass  filtering  with  depth.  For  temporal  frequencies  that  are  high  enough  compared  to  the 
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spatial  frequency,  i.e.  on  coarse  spatial  scale,  the  exponent  becomes  imaginary  and  therefore 
corresponds  to  a  time-delay  caused  by  the  wavefront  traveling  at  a  finite  speed,  (cp  or  csh). 
For  finer  spatial  details,  the  exponent  serves  the  same  purpose  as  in  the  static  case,  i.e.  the 
increased  spatial  low-pass  filtering  (blurring)  of  signals  with  depth. 

Neglecting  inertia:  For  further  analysis  of  the  transfer  function,  it  is  helpful  to  renormalize 
and  measure  distances  and  displacements  in  units  of  shear- wave  wavelengths.  Let  s  =  ju, 
X sh  =  cshjuj  and  hence  z  =  z\sh ,  ujx  =  Cjxj\sh.  The  can  then  be  canceled  out  from  the 
transfer  function  which  then  becomes 

nt  [{(%  -  l/2)e~‘**  -  age-"**] 

Uz  2/z  [(w2  -  1/2)2  -  w2nin2] 

with  h\  =  (<Z'l  -  /J-2)1/2,  and  n2  =  (u>2  -  l)1/2.  We  observe  that  there  are  no  singularities 
and  all  the  terms  are  real  in  the  transfer  function  for  Cjx  >  1.  This  can  be  interpreted  as 
spatial  frequencies  with  Cjx  >  1  are  not  involved  in  energy  transfer  away  from  the  source  and 
that  the  halfspace  behaves  simply  as  a  spring  for  those  frequencies.  The  denominator  however 
vanishes  for  Cjx  =  cnjcsh ,  since  it  is  identical  to  Eq.  (2.38)  used  to  obtain  cn/cSh .  That  the 
halfspace  behaves  as  a  spring  for  some  frequency  range  becomes  important  when  one  considers 
the  conditions  under  which  inertia  can  be  neglected,  especially  in  cases  where  the  solution  in  Eq. 
(2.41)  is  applied  to  finite  domains  as  will  eventually  happen  in  our  case.  Let  x  €  [— L/2,  Lj 2], 
then  only  a  discrete  set  of  spatial  frequencies  will  give  valid  solutions  as  discussed  in  section 
2.4.  The  lowest  nonzero  frequency  is  of  particular  interest  to  us  and  the  modeshapes  (assuming 
simple  boundary  conditions  at  x  —  ±L/ 2)  corresponding  to  that  have  been  plotted  in  Fig.  2-6. 
We  see  that  the  longest  possible  wavelength  is  Amax  =  2 L  which  is  equivalent  to  min(o;x)  =  7 x/L. 
We  therefore  have  that  Cjx  <  cn/csh  is  satisfied  if 


giving  us  a  condition  as  to  when  we  can  ignore  the  effects  of  inertia  for  the  semi- infinite  halfspace. 


4.5  Viscoelasticity 

Theory  of  viscoelasticity  extends  the  theory  of  elasticity  to  include  cases  where  the  stress-strain 
relation  behavior  is  rate  dependent.  It  all  boils  down  to  the  choice  of  timescale,  and  for  shorter 
intervals  of  interest,  it  may  be  necessary  to  account  the  rate  dependent  behavior  of  materials. 
The  description  linear  viscoelastic  material  applies  to  materials  for  which  the  rate  dependent 
behavior  can  be  adequately  described  by  expressing  the  stress-strain  relation  in  the  viscoelastic 
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Figure  2-6:  Shows  the  modeshape  for  the  lowest  eigenfrequeneies  of  a  domain  of  length  L.  The  maximum 
wavelength  is  Amax  =  2 L  and  hence  the  lowest  non- zero  spatial  eigenfrequency  is  min(wx.)  =  tt/L. 


form  of  Hooke’s  law, 


<r(s)  =  Q(s)e(s). 


Where  a(s)  and  e(s)  are  the  Laplace  transforms  of  the  stress  and  the  strain,  respectively.  Q(s) 
is  called  the  Operational  relaxance  and  alternatively  Q(s)/s  is  called  the  Operational  impedance. 
The  main  advantage  in  using  the  viscoelastic  form  of  Hooke’s  law  is  that  we  can  then  apply  the 
correspondence  principle.  According  to  the  principle,  if  an  elastic  solution  to  a  stress  analysis 
problem  is  known,  substitution  of  the  appropriate  Laplace  transforms  for  the  material  constants 
employed  in  the  elastic  solution  furnishes  the  viscoelastic  solutions  [83]. 


Hence,  the  effects  of  viscoelasticity  can  be  incorporated  into  the  elastodynamic  solution 
in  (2.41)  simply  by  allowing  the  material  constants  to  be  frequency  dependent,  i.e.  //  and  (3 
become  p(s)  and  3(s).  Similarly  the  effects  of  viscoelasticity  can  included  in  the  static  solution 
in  Eq.  (2.30)  or  in  Eqs.  (2.32)  and  (2.33). 


The  material  constants  are  in  general  both  independently  rate  independent,  p  describing 
shear  motion  and  3  the  dilational  behavior.  It  is  however  often  assumed  that  only  //.  is  rate 
dependent.  This  assumption  is  reasonable  for  materials  that  are  near  being  incompressible, 
(u  >  0.4)  which  is  the  case  for  polymers  and  other  materials  we  intent  to  apply  this  theory  to. 
The  simplest  model  that  captures  the  behavior  of  a  linearly  viscoelastic  solid  is  the  3-parameter 
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(a)  3-Parameter  Standard  Linear  Solid 


(b)  Response  of  stress  to  a  step  in  strain 


Figure  2-7:  (a)  Shows  a  diagram  of  the  3-parameter  linear  standard  solid  model  using  spring  and  damper 
elements,  (b)  Shows  the  response  in  stress  to  a  step  change  in  strain  for  the  3-parameter  linear  solid 
defined  in  Fig.  2-7.  The  relation  between  the  parameters  is  that  fi0  =  +  ^2)1  Q'  =  (^’1  +  ^2)/T'2, 

and  tr  —  b\/(ki  k-2). 


standard  linear  solid.9  For  this  model  the  transfer  function  for  /a  is, 


H(s)  -  ^ 


1  4-  trs 
1  +  ~trs 


where  no  is  the  steady  state  response,  tr  the  retardation  time  and  a  is  the  overshoot  of  a  step 
response.  A  schematic  representation  of  the  model  as  dampers  and  springs  is  shown  in  Fig. 
2-7a,  and  the  step  response  in  Fig.  2-7b,  i.e.  the  response  in  stress  to  a  step  in  strain. 


This  model  can  be  generalized  to  model  materials  with  more  than  one  time  constant.  The 
transfer  function  for  the  generalized  model  is 


m(s) = 40  n 


1  +  tr-  S 


.  ,  1  +  —tr-S 

?  =  1  Qi  '< 


The  transfer  function  approach  to  viscoelasticity  provides  a  powerful  tool  for  the  modeling 
of  the  viscoelastic  properties  of  different  materials.  One  still  needs,  however,  to  estimate  the 
parameters  in  the  transfer  function  for  each  material  and  under  conditions  that  are  as  close  as 
possible  to  the  conditions  under  which  the  model  is  to  be  applied. 


9A1so  called  Kelvin  model. 
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5.  Conclusions 


In  this  chapter  we  have  given  the  basic  theory  of  integral  transforms  and  shown  how  they  can  be 
applied  to  problems  in  contact  mechanics.  We  have  shown  how  the  solutions  to  contact  problems 
involving  various  structures  can  be  presented  as  transfer  functions.  The  transfer  function  takes 
as  inputs  the  boundary  conditions  but  gives  the  desired  solution  as  output.  Analysis  of  the 
transfer  function  has  also  been  shown  to  give  us  the  qualitative  properties  and  behavior  of  the 
solution,  irrespective  of  the  exact  from  of  the  input.  The  transfer  function  approach  has  further 
been  shown  to  be  equally  applicable  to  problems  involving  static,  dynamic  and  viscoelastic 
effects,  and  further  analysis  of  the  solutions  has  given  us  conditions  under  which  inertial  effects 
can  be  ignored  and  the  quasi-static  solutions  are  valid. 


41 


42 


Chapter  3 


Static  Identification 


1.  Introduction 


Manual  exploration  and  manipulation  of  objects  in  an  environment  is  important  to  both  hu¬ 
mans  and  robots.  Human  performance  of  these  tasks  is  enabled  by  the  haptic  system  consisting 
of  tactile  and  kinesthetic  sensory  systems  together  with  a  motor  control  system.  Similar  classi¬ 
fication  is  applicable  to  the  sensory  and  motor  systems  of  robots.  A  detailed  and  quantitative 
understanding  of  the  underlying  dynamics,  information  flow,  and  control  strategies  will  benefit 
investigations  of  both  human  haptics  and  the  development  of  machine  haptics.  It  is  especially 
valuable  in  the  development  of  haptic  interfaces  through  which  humans  can  interact  manually 
with  teleoperated  systems  or  computer  generated  virtual  environments  [74].  Although  the  prin¬ 
ciples  of  operation  of  man-made  devices  are  quite  different  from  those  of  humans,  the  constraints 
on  the  performance  of  these  haptic  tasks,  such  as  the  laws  of  physics  governing  the  mechanics  of 
contact  and  the  presence  of  friction  and  gravity  are  the  same  for  both.  In  addition,  the  type  of 
tactual  sensory  information,  their  processing  and  the  computation  of  the  required  control  action 
are  sufficiently  similar  for  the  two  systems  that  the  common  aspects  of  information  processing 
can  be  functionally  separated  from  the  hardware  implementations.  Therefore  a  computational 
theory  of  haptics  that  investigates  what  kind  of  information  is  necessary,  and  how  it  has  to  be 
processed  in  order  to  successfully  complete  a  desired  haptic  task  can  be  common  to  humans 
and  robotic  systems. 

It  has  been  recognized  for  some  time  that  robots  with  dextrous  hands  need  ‘tactile  sense’ 
in  order  to  successfully  explore  and  manipulate  objects  in  their  environment.  Reviews  of  the 
various  tactile  sensors  designs  have  been  given  by  several  authors,  (see  for  example  [59,  68,  32].) 
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Because  of  the  natural  constraints  imposed  by  this  mode  of  sensing,  the  engineering  designs 
have  had  to  follow  nature  in  overall  configuration.  By  definition,  tactile  sensing  is  achieved 
through  direct  contact  with  objects,  and  therefore  a  ‘skin’  is  necessary  to  protect  the  sensors 
from  physical  damage.  The  requirements  that  the  skin  should  be  soft  comes  from  the  need  to 
have  (1)  regions  of  contact  within  which  skin  surface  conforms  to  the  object  surface  (instead  of 
point  or  line  contact  that  occurs  between  two  rigid  objects)  (2)  significant  deformation  within 
the  medium  so  that  the  sensors  are  activated  and  have  enough  resolution.  If  the  substrate 
material  on  which  the  sensors  and  the  skin  rest  is  also  soft,  then,  in  addition  to  the  above, 
better  prehension  stability  can  be  achieved.  The  fact  that  a  point  contact  has  no  torsional 
resistance  and  that  its  stability  is  highly  sensitive  to  local  aberrations  implies  that  compliant 
skin  gives  better  prehension.  Thus,  although  robotic  tactile  sensors  themselves  might  differ  in 
their  operation,  depending  upon  whether  they  respond  to  changes  in  conductance,  capacitance, 
contact  area,  light  intensity,  etc.,  within  each  sensing  cell  caused  be  local  mechanical  distortions, 
the  overall  configuration  of  all  the  designs  is  that  of  mechano-sensitive  transducers  embedded 
in  a  deformable  medium. 

Two  stages  in  tactual  information  processing,  encoding  and  decoding  make  up  the  identifica¬ 
tion  problem ,  which  is  so  called  because  it  involves  the  determination  of  the  state  of  the  object 
and  its  relationship  to  the  tactual  sensory  system.  The  system  discussed  here  differs  in  some 
ways  from  traditional  dynamic  systems.  The  most  prominent  difference  is  that  it  has  both  spa¬ 
tial  and  temporal  variables,  whereas  traditional  lumped  parameter  systems  only  have  temporal 
ones.  In  both  natural  and  man-made  tactile  systems  the  encoding  process  is  the  transduction  of 
the  mechanical  stimulus  into  electrical  signals  that  contain,  in  a  coded  form,  information  about 
the  stimulus  features.  The  inverse  problem,  i.e.  the  decoding  process  or  the  calculation  of  the 
surface  shape  or  surface  tractions  from  the  transduced  signals  is  in  many  respects  significantly 
more  difficult.  The  encoding  is  essentially  a  low-pass  filtering  process  which  indicates  that  the 
decoding  will  be  ill-posecl  in  the  sense  that  the  higher  frequency  components  of  the  signals  will 
be  prone  to  noise  contamination  and  special  methods  have  to  be  applied  for  successful  decoding. 

Several  papers  can  be  found  in  the  robotics  literature  on  the  mechanistic  modeling,  encoding 
and  decoding  of  surface  conditions  from  subsurface  measurements  using  tactile  sensors.  The 
earlier  papers  dealt  primarily  with  the  mechanistic  model  and  only  the  encoding.  In  [35]  the 
problem  is  approached  from  a  physiological  point  of  view  and  in  [20]  from  robotic  point  of 
view.  Both  papers  use  the  plane  strain  assumption  and  give  solutions  to  the  2-dimensional 
(2D)  encoding  problem.  The  decoding  problem  is  first  addressed  in  [73,  63],  both  of  which 
solve  only  the  2D  plane  strain  case.  [73]  uses  a  frequency  domain  approach  whereas  a  neural 
network  approach  is  used  in  [63].  The  plane  strain  assumption  is  dropped  in  [9]  but  the  analysis 
is  limited  to  the  properties  of  a  specific  photoelastic  sensor.  In  [18]  again  plane  stress  or  plane 
strain  is  assumed,  the  effects  of  finite  thickness  and  the  stability  of  the  grasp  are  investigated. 
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The  mathematics  of  the  decoding  problem  are  investigated  using  matrix  regularizing  operators 
in  [65],  where  the  loading  is  assumed  to  be  axisymmetric.  The  combined  results  of  the  above 
papers  are  detailed  in  [19]  which  is  the  most  thorough  treatment  so  far  on  the  infinite  half¬ 
space  and  analytical  solutions  of  the  encoding  and  the  decoding  problem.  [19]  further  gives  the 
3-dimensional  (3D)  half-space  solution  only  to  a  normal  point  load  in  appendix.  Additionally, 
finite  element  analysis  is  used  in  [72],  where  the  effects  of  shear  loads  on  sub-surface  strain  and 
stress  are  demonstrated  graphically,  and  also  in  [17]  which  investigates  the  differentiation  of 
shapes  and  the  inherent  difficulties  that  arise  because  of  the  blurring  effect  of  the  skin.  In  [28] 
the  encoding  and  decoding  of  a  stress  rate  sensor  are  investigated  and  a  scalar  Wiener  inverse  is 
suggested  for  the  decoding.  Finally,  in  [60]  the  2D  plane  strain  solution  for  a  cylindrical  finger 
is  given. 

In  summary,  none  of  the  papers  above  treats  either  the  encoding  or  the  decoding  problem 
fully.  All  the  results  published  until  now  on  the  reconstruction  of  general  surface  loads  and  pro¬ 
files  from  subsurface  tactile  signals  do  not  solve  the  full  3D  problem.  The  effects  and  interaction 
of  shear  loads  and  the  selection  of  appropriate  sensor  signals  has  not  been  investigated  before. 
In  addition,  neither  the  decoding  problem  has  been  fully  addressed  nor  has  the  multivariate 
Wiener  inverse  been  applied.  In  this  chapter  we  will  address  mainly  the  mechanistic  modeling, 
encoding,  and  the  decoding  problems.  The  contributions  of  this  chapter  can  be  itemized  as 
follows: 


1.  A  complete  solution  to  the  3D  half-space  problem  is  given  for  arbitrary  3D  static  loading 
conditions  and  any  linearly  elastic  material,  instead  of  the  2D  solutions  used  before.  The 
solution  also  applies  to  pseudo  dynamic  problems  which  can  be  described  as  a  sequence 
of  elasto-static  problems  (such  as  the  estimation  of  the  onset  of  slip). 

2.  The  solution  is  formulated  using  transfer  function  matrices  in  the  domain  of  spatial  fre¬ 
quency  which  allows  us  to  look  at  the  qualitative  properties  and  general  behavior  of  the 
solution,  irrespective  of  particular  cases  of  the  loading  or  object  shape.  It  also  enables 
efficient  numerical  implementation  based  on  the  theory  of  linear  signal  processing. 

3.  Symmetry  and  signal  bandwidth  arguments  are  invoked  to  select  the  optimal  transducing 
signals  or  types  of  sensors.  For  decoding,  we  can  then  apply  the  signal  processing  and 
regularization  theory  in  a  manner  similar  to  what,  has  been  done  in  the  development  of 
computational  theory  of  vision. 

4.  Numerically  stable  multivariate  regularizing  solution  methods  for  the  decoding  problem 
are  delineated.  The  ill-posedness  of  the  decoding  process  is  analyzed  and  it  is  shown  as 
to  how  the  extensively  developed  methods  of  image  restoration  can  be  applied  to  avoid 
the  problems  associated  with  the  ill-posedness. 
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In  Section  2,  we  formulate  and  solve  the  3D  half-space  elastic  model,  which  serves  as  a 
description  of  the  encoding  process.  For  the  sake  of  clarity,  the  mathematical  expressions  are 
detailed  in  Appendix  A.  In  Section  3,  we  analyze  the  decoding  problem,  first  as  an  ideal 
problem  without  any  real  world  effects  such  as  noise,  following  which  we  incorporate  the  effects 
of  sampling  and  random  additive  noise  into  the  solution.  This  constitutes  an  algorithm  for  the 
decoding  process.  In  Section  4,  we  give  example  numerical  simulations  of  the  methods  proposed 
in  Section  3.  Conclusions  are  given  in  section  5. 


2.  The  Encoding  Problem 

From  a  theoretical  viewpoint,  tactile  sensor  response  to  any  given  stimulus  can  be  predicted  if 
we  have  the  following: 

1.  A  mechanistic  model  of  the  deformable  medium  which  can  be  used  to  calculate  reliably 
the  stresses  and  strains  at  each  point  in  the  medium  for  a  given  mechanical  stimulus. 

2.  A  model  of  the  sensor  that  provides  the  relationship  between  the  relevant  stimulus,  i.e.,  a 
particular  combination  of  stresses  and  strains  in  the  local  neighborhood  of  a  sensor  that 
it  is  responsive  to  and  the  sensor  response. 

For  example,  an  appropriate  mechanistic  model  of  the  human  skin  and  the  soft  tissues  under¬ 
neath  enables  us  to  investigate  one  of  the  important  questions  in  cutaneous  neurophysiology, 
that  of  identifying  the  relevant  stimulus  that  causes  each  receptor  type  to  respond  [35,  75] 


2.1  Mechanistic  Analysis 

In  the  robotic  tactile  sensing  literature,  contact  analysis  problems  are  frequently  simplified  by 
posing  them  as  2D  plane  strain  problems  in  the  elastic  half-space  [70,  23,  36].  The  surface  load 
or  profile  is  then  assumed  not  to  change  along  one  dimension  on  the  planar  surface,  for  example, 
as  in  the  case  of  a  line  load.  Almost  all  the  results  in  the  literature  on  the  analysis  of  tactile 
sensing  have  made  this  assumption.  Clearly,  this  is  a  serious  limitation.  In  this  section  we 
will  give  the  more  general  solution  when  3D  loads  are  arbitrarily  distributed  over  the  contact 
region. 

The  assumption  that  the  problem  can  be  posed  as  a  contact  problem  in  the  semi-infinite 
half-space  still  remains,  mainly  for  the  sake  of  analytical  tractability.  We  can,  however,  partially 


justify  that  assumption  based  on  numerical  solutions  obtained  for  finite  models  [13,  75]  which 
show  that  the  effects  of  finite  extent  are  minimal  when  the  contact  region  is  small  relative  to 
the  surface  area  of  the  medium.  The  semi-infinite  half-space  solution  is  a  good  approximation 
in  such  cases  because  the  stress/strain  response  to  a  load  decreases  exponentially  with  distance 
from  the  load. 

The  problem  is  set  up  so  that  the  boundary  surface  coincides  with  the  rcy-plane  and  the 
positive  2-axis  points  into  the  medium  (See  Fig.  3-1).  The  material  is  assumed  to  be  lin- 


Figure  3-1:  Shows  the  semi-infinite  halfspace  with  surface  tractions  and  subsurface  sensors. 
The  x-  and  y-dimensions  extend  to  infinity  in  both  negative  and  positive  directions.  The  z- 
dimension  extends  from  0  to  positive  infinity.  The  tractions  are  distributed  over  the  contact 
area,  each  point  being  subjected  to  3-dimensional  loading  f(x,y)  =  [fx,  fy,  fz}T. 

early  elastic  and  homogeneous.  It  is  also  assumed  to  extend  to  infinity  in  both  negative  and 
positive  x ,y  directions  and  in  the  positive  2  direction.  The  surface  load  is  described  by  the 
traction  vector  field  f{x,y )  =  [fx(x,y),  fy(x,y ),  fz{x,  y)]T  and  the  surface  displacement  field 
as  uo(x,y)  =  [ux(.r,  y,  0),  uy(x,y. 0),  uz(x,  y,  0)]’r.  In  our  formulation,  the  boundary  conditions 
can  be  forces,  displacements  or  mixed,  i.e.  when  any  3  elements  from  both  /  and  uq  are  given. 
We  intend  to  present  the  solution  in  the  space  of  spatial  frequencies  and  use  transfer  functions 
to  simplify  solutions  and  gain  further  insight. 


2.2  Transfer  Function  Matrices  (TFM) 

Many  equivalent  routes  to  the  solution  of  the  encoding  problem  exist,  but  the  simplest  and 
most  direct  way  for  our  purposes  is  to  transform  the  fundamental  differential  equations  of 
equilibrium  and  the  stress-strain  relationship  equations  (Generalized  Hooke’s  law)  using  the 
Fourier  transform  over  (x,y)  space.  Examples  of  this  approach  can  be  seen  in  [16]  and  [71]. 
The  derivation  of  the  transfer  function  matrix  form  of  the  solution  is  given  in  Appendix  A. 

Similar  to  what  is  done  in  control  theory,  we  employ  here  the  concept  of  a  transfer  function 
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matrix  (TFM).  The  transfer  function  matrices  we  give  here  have  the  Fourier  transform  of 
the  components  of  either  f(x,y)  or  uq (x,y)  as  the  inputs  and  any  displacements,  strains  or 
stresses  in  the  medium  as  the  outputs.  Each  entry  of  the  TFM  would  then  correspond  to  the 
contribution  of  a  particular  input  component  to  the  corresponding  output  component.  For 
example  the  entry  in  the  TFM  for  the  ^-displacement  resulting  from  the  x-load,  or  uz(x,y,z) 
we  would  write 

Hz(u)x,UJyi  -  TUX  Wy,  ^y) 

to  obtain  the  Fourier  transform  uxz (u>x,uyi  z)  of  uz(x,y,z)  from  the  transfer  function 
Tu* f(ojx,u)yi  z)  and  the  Fourier  transform  fx(uix,ujy)  of  the  surface  shear  load  fx(x,y).  Here, 
Lox  and  iOy  correspond  to  spatial  frequencies  in  the  x  and  y  directions  respectively.  Then  we 
can  write 

^  =  Tuff  e  —  Teff  (7  =  Taff  pn  —  TPnff 

(3.1) 

^  —  TUUq11q  £  —  TcUoUq  a  —  T<yUo1lo  pn  —  TpnU0U() 


where  the  subscript  /  indicates  that  surface  tractions  are  inputs  and  the  subscript  uq  indicates 
that  surface  displacements  are  inputs.  As  an  example,  we  give  here  Tsuoi  i-e-  the  case  when 
the  surface  displacements  («o)  are  the  inputs  and  the  subsurface  strains  (e)  are  outputs.  The 
other  TFMs  are  given  in  Appendix  A. 


e"*n 

e"°  “  (3  -  4i/)Q 


jux  ((3  -  4i/)n  -u %z) 

— jujxujyZ 

—jQujx  (1  —  flz) 


'y  ((3  -  4 v)n  -  2 J*xz) 


2W 

-i((3-4i/)fr2  +  w2(i-2ns)) 

-^Uxujy{l  -  2 zQ) 

,1/2 


-JU&yZ 

juy  ((3  -  4v)Sl  -  w2z) 

—  jQ.UJy  (1  —  Oz) 

((3  -  4i/)n  -  2 x2z) 

-^UJxU!y(  1  -  2zfl) 


fix2 


l 'iXyZ 


-n2(2(l  -  2is)  +  ftz) 


fixx  Xy  .. 


jttux(l  —  2v)  +  Qz) 


((3  -  4i/)n2  +  x2(l  -  2f lz))  jflux{  1  -  2v)  +  Hz) 


(3.2) 

.'here  fl  =  (x2  +  x2)1^,  y  is  the  shear  modulus,  and  v  is  the  Poisson  ratio.  TFMs  for  mixed  cases 
can  also  be  easily  constructed  by  choosing  appropriate  elements  from  the  TFMs  above. 


The  equations  in  (3.1)  represent  the  solution  for  the  encoding  problem,  because  given 
/(x,  y),  il0(x ,  y)  or  a  combination  of  3  elements  from  the  two,  we  can  predict  u(x,  y ,  z),  e(x,  y,  z), 
a(x,y,z)  and  pn(x,y,z).  If  the  relevant  stimulus,  i.e.  a  particular  combination  of  stresses  and 
strains  in  the  local  neighborhood  of  the  sensor  that  it  is  responsive  to,  is  known,  we  can  predict 
the  sensor  response.  For  example,  if  the  sensor  response  is  linearly  related  to  the  strain  £zz, 
then  it  is  the  relevant  stimulus  for  that  sensor.  On  the  other  hand,  if  the  relevant  stimulus 
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needs  to  be  identified,  we  can  use  Eq.  (3.1)  to  generate  hypotheses  that  can  be  tested  against 
the  empirically  measured  sensor  response,  as  is  done  in  [35,  75]. 


3.  The  Decoding  Problem 


Tactile  information  is  obtained  with  mechanosensors  embedded  in  a  deformable  medium  that 
is  in  contact  with  a  shaped  object.  Mechanosensors,  embedded  within  the  medium  can,  at 
best,  provide  information  on  skin  surface  shape  and  surface  traction  distribution.  From  this 
information,  contact  region,  shape  of  the  object  and  the  contact  state  need  to  be  inferred. 
The  problem  at  hand  is  therefore  the  decoding,  i.e.  reconstruction  of  the  surface  shape  of 
the  medium,  the  tractions  on  the  surface  as  well  as  the  region  of  contact  from  sub-surface 
mechanosensory  information. 

The  approach  we  have  chosen  is  to  represent  the  solution  in  the  spatial  frequency  domain,  in 
terms  of  transfer  functions  as  in  Eq.  (3.1).  Equivalently,  using  superposition  of  the  point  load 
solution,  the  solutions  can  also  be  presented  as  convolution  integrals.  The  transfer  function 
approach  gives,  however,  solutions  whose  properties  are  more  transparent  and  some  important 
statements  can  be  made  regarding  their  qualitative  properties  (such  as  the  symmetry  and 
bandwidth  properties  of  subsurface  stress  or  strain  components),  irrespective  of  the  exact  form 
of  the  loading.  This  approach  also  enables  the  use  of  efficient  tools  of  linear  signal  processing, 
such  as  the  Fast  Fourier  Transform. 

This  section  is  organized  as  follows.  In  section  3.1  we  discuss  the  constraints  that  the  task 
of  tactile  identification  puts  on  the  selection  and  the  number  of  necessary  sensor  signals,  and 
furthermore  identify  possible  candidates  of  sensor  signal  combinations.  In  section  3.2  we  use 
symmetry  arguments  to  analyze  the  TFMs  of  the  previously  identified  candidates  in  order  to 
find  an  optimal  combination.  Finally,  in  section  3.3  we  show  how  an  optimal  multivariate 
Wiener  inverse  can  be  used  to  reconstruct  the  surface  signals  from  measurements  of  the  sensor 
signal  in  the  presence  of  noise  and  sampling  errors. 
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3.1  Sensory  Information  Necessary  for  Reconstruction 


The  most  obvious  signals  to  reconstruct  are  the  deformed  shape  of  the  surface  of  the  medium  and 
the  spatial  distribution  of  surface  tractions.  The  extents  of  the  contact  area  need  to  be  assessed 
but  the  surface  displacements  are  necessarily  continuous  at  the  edge  of  the  contact  area,  so  they 
provide  limited  information.  Discontinuities  in  slope  or  curvature  of  surface  shape  however,  do 
not  reliably  indicate  contact  borders  as  an  object  surface  within  the  contact  region  could  have 
similar  discontinuities  and  therefore  be  mistakenly  identified  as  borders  of  the  contact  region. 
Each  of  the  regions  where  the  surface  tractions  (e.g.  normal  pressure)  are  non-zero  indicates 
the  region  of  contact  with  the  object.  Within  these  regions  the  object  shape  is  the  same  as 
the  shape  of  the  deformed  medium  surface.  In  addition,  the  surface  traction  distribution  can 
directly  give  valuable  information  on  the  stability  of  grasp  and  the  type  of  contact,  be  it  static 
or  vibrating. 

Some  interesting  properties  of  the  TFM  can  be  obtained  when  dynamic  systems  theory  is 
applied  to  them.  In  what  follows  we  assume  that  the  sensors  pick  up  a  strain  component  or 
mean  normal  stress  and  further  that  the  components  of  the  strain  are  expressed  in  the  global 
coordinate  system  ( x,y,z ).  If  the  medium  is  incompressible,  uniform  normal  pressure  on  its 
surface  does  not  cause  any  strains  within.  Therefore,  measurements  of  strains  only  are  not 
sufficient  to  fully  reconstruct  surface  force  distributions,  as  the  mean  normal  stress  component 
will  not  be  detectable  in  the  strains  for  incompressible  materials.  This  is  the  main  reason  for 
the  use  of  the  mean  normal  stress  as  a  variable  that  needs  to  be  sensed. 

If  no  shear  loads  are  present,  or  if  the  relationship  between  normal  and  shear  loads  is  known, 
only  one  component  of  strain  needs  to  be  sensed.  However,  if  shear  independent  of  the  normal 
load  is  present,  the  more  general  TFM  formulation  must  be  employed,  as  the  surface  tractions 
must  be  obtained  by  the  inversion  of  some  of  the  equations  in  Eq.  (A. 7)  and  (A. 9)  in  Appendix 
A.  What  and  how  many  components  of  the  strain  tensor  are  needed  depend  on  the  assumptions 
we  make  on  the  loading  conditions.  This  can  be  detailed  in  the  following  cases. 


1.  No  shear  loading:  In  this  case  the  normal  pressure  distribution  and  surface  displace¬ 
ments  need  to  be  reconstructed.  Sensing  of  any  strain  or  the  mean  normal  pressure  is 
generally  sufficient,  but  the  need  for  rotational  symmetry  about  the  2-axis  makes  either  of 
Pn  or  ezz  the  natural  choice.  Using  pn  will  enable  the  calculation  of  the  uniform  pressure 
component  of  the  total  load  which  no  combination  of  strains  can  provide  if  u  =  0.5  (i.e. 
incompressible  material).  However,  the  combination  of  pn  and  ezz  can  also  be  used  to 
obtain  increased  spatial  frequency  bandwidth  (see  section  3.2),  in  addition  to  obtaining 
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the  full  load  distribution. 


2.  Relationship  between  normal  loads  and  shear  loads  is  fully  or  partially  known: 
This  is  the  case  for  example,  if  the  object  is  sliding  over  in  an  unknown  direction  over 
the  surface  under  Coulomb  friction  conditions,  or  if  the  relationship  between  normal  and 
shear  loads  is  unknown  but  the  direction  of  the  shear  load  is  known.  Then  any  two 
subsurface  stress  and  strain  components  need  to  be  sensed  to  fully  reconstruct  all  three 
surface  traction  or  displacement  components,  and  a  natural  choice  would  again  be  pn  and 
£zz  as  above. 

3.  All  tractions  unknown:  All  3  traction  components  and  surface  displacements  need 
to  be  reconstructed.  Considering  only  combinations  symmetric  in  (x,  y)  the  following 
combinations  become  candidates:  (a)  (sxy,exx,s:yy),  (b)  (£xx,£yy,£zz),  (c)  (£zz,  £xz,  £yz), 

(d)  (f  xy ,  4xz  j  &yz ) )  (®)  {Pni  ^zzi  ^xy  )»  (f )  Pm^xx^yyi  3-Ild  (g)  ijPni^xzi^yz')- 


We  will  shortly  show  how  we  can  make  an  optimal  choice  among  these  component  sets. 

Case  3  is  the  most  general  one  and  it  indicates  that  any  sensor  population  that  is  suit¬ 
able  for  reconstructing  arbitrary  loading  conditions  must  measure  at  least  three  independent 
stress /strain  components  (out  of  a  total  of  six).  More  components  could  be  used  to  get  re¬ 
dundant  measurements,  and  hence  more  reliable  estimates,  but  here  we  will  only  look  at  the 
minimal  set. 


3.2  Analysis  of  TFM  Properties 


We  now  try  to  use  the  properties  of  the  TFM  to  identify  the  sensor  combination  that  is  most 
suitable  for  tactile  sensing.  Particulary  we  will  a  look  at  the  rank  of  the  TFM  and  its  bandwidth. 


Rank  of  the  TFM 


Since  we  intend  to  invert  the  TFM,  we  need  to  consider  its  rank,  which  is  a  key  property  and 
has  to  be  full  for  an  inversion  to  be  possible.  If  the  TFM  loses  rank,  information  is  lost  as  the 
output  will  span  fewer  dimensions  than  the  input.  The  conditions  under  which  the  TFM  looses 
rank  are  closely  related  to  zero-properties  of  dynamic  systems,  just  as  the  conditions  under 
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which  the  TFM  becomes  infinite  are  related  to  the  poles.  The  TFMs  in  our  case  are  however 
always  finite  and  therefore  do  not  possess  any  pole-like  properties,  which  is  to  be  expected  from 
a  linear  mechanical  system  in  stable  equilibrium. 

The  TFM  looses  rank  when  its  determinant  becomes  zero.  When  the  traction  vector  /  is 
the  input,  the  determinants  of  the  TFM  candidates  detailed  in  section  3.1  are 


det(Tfxj/i£ixi£!/S() 

=  0 

(a) 

det {Texx,£yy,£zz  ) 

=  -(l-2V)^e-n 

(b) 

det  (T£zZ:SjCZt£yz ) 

=  — (1  —  2u)e~z^ 

(c) 

det(T£ij)i£xzi£:!/J. ) 

(d) 

det  {TpnteZZ}€iy ) 

2(1  +  v)  -  u2y  _zQ 

3  n2  6 

(e) 

det(T;  pn,exx,£yy) 

2(1  +  f)  2 UXUy  _.Q 

3  Q2  6 

(f) 

det(TPn  ,eIZl£yz) 

2(1  +  v)  n 

3 

(g) 

We  see  that  the  first  one  is  always  rank  deficient.  If  the  medium  is  incompressible,  i.e.  if 
1  —  2v  =  0,  then  the  next  three  TFMs  loose  rank.  As  discussed  in  section  3.1  the  mean 
normal  stress,  pn.  must  be  among  the  sensor  signals  in  the  case  of  an  incompressible  material. 
Otherwise  the  component  of  the  surface  load  that  contributes  to  the  mean  normal  stress  cannot 
be  reconstructed.  A  special  case  is  a  solid  subjected  to  uniform  pressure.  No  strains  will  be 
detected  inside  the  solid,  although  the  pressure  can  be  arbitrary. 

We  also  note  that  the  TFM  looses  rank  for  candidate  (e)  when  a )x  =  ±a iy  and  for  candidate 
(f)  when  ux  =  0  or  ioy  =  0.  That  the  TFM  looses  rank  for  certain  directions  in  the  (ljx,  ujy) 
plane,  means  that  the  inputs  cannot  be  reconstructed  along  those  same  lines.  In  the  case  of 
candidate  (f)  this  means  that  the  cumulative  load  distribution  in  either  x  or  y  direction  cannot 
be  reconstructed.  For  example,  if  u>y  =  0,  in  frequency  space  we  have  f(ux,0)  =  g(ux)  and 
g(x)  =  f(x,y)dy,  which  is  the  cumulative  load  distribution  along  x.  Therefore,  since  the 
TFM  for  candidate  (f),  T(f),  is  rank  deficient  for  uy  —  0,  the  component  containing  information 
on  the  cumulative  distribution  of  the  input  along  x  will  be  lost  in  the  encoding.  The  same 
arguments  apply  to  candidate  (e),  with  the  coordinate  system  rotated  by  ±45°  around  the 
2-axis.  Candidate  (g)  is  the  only  one  that  does  not  loose  rank  and  therefore  has  the  desirable 
combination  of  p„  and  strains  that  the  sensors  need  to  transduce  for  decoding. 
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Remark:  The  comments  above  are  also  valid  when  the  displacement  vector  uo  or  any  combi¬ 
nation  of  tractions  and  displacements  are  the  inputs,  since  TUUo  =  TufTfuo  and  the  TFM  that 
relates  u o  and  /  (Eq.  (A. 10)  in  the  Appendix)  has  the  determinant, 


det  Tuof  = 


2(1  +  3  1 

3  —  4i/  fl3 


(3.3) 


which  does  not  alter  the  rank  properties  of  each  of  the  TFMs  of  candidates  (a)-(g)  above. 


Bandwidth  of  Spatial  Frequency 


A  frequently  used  property  of  dynamic  systems  is  the  bandwidth  of  the  components  of  the 
TFM  which  measures  how  the  attenuation  of  each  signal  in  the  medium  depends  on  frequency. 
High  bandwidth  would  mean  that  reconstruction  is  possible  with  finer  details  and  is  therefore 
desirable.  However,  since  the  frequency  is  2D,  the  bandwidth  will  also  depend  on  orientation  in 
the  (a;,  y)  plane,  which  in  turn  can  be  used  to  exclude  sensor  combinations  which  lack  rotational 
symmetry  about  the  2-axis. 


To  compare  the  bandwidth  of  the  various  combinations  (a)-(g)  listed  above,  we  look  at  the 
power  spectrum  of  the  TFMs,  a  line  at  a  time.  That  is,  the  power  spectrum  of  each  sensor 
signal  is  the  norm  of  the  corresponding  line  vector  in  the  TFM.  For  example 


V£ 


xz 


TF- 


Tel, 


but  this  implicitly  assumes  that  the  all  the  traction  imputs  are  weighted  equally.  Weighing  the 
normal  traction  differently  than  the  shear  traction  does  not  alter  the  result  of  the  following 
arguments,  as  long  as  the  shear  tractions  are  weighted  equally.  The  power  spectrum  for  the 
possible  sensory  signals  can  be  seen  in  Figs.  3-2-3-7.  The  figures  show  the  power  spectrum  as 
density  plots  with  darker  and  lighter  areas  indicating  frequencies  where  attenuation  is  low  and 
high,  respectively.  Also  shown  along  the  border  of  the  density  plot,  are  projected  profiles  of  the 
power  spectrum  along  the  lines  ux  =  0,  uy  —  0,  and  ux  =  uy.  Figures  3- 2-3-7  show  the  power 
spectrum  respectively  for  pn,  £xx,£zz,£xy,exz  and  candidate  (g)  (VPn  +V£xz  +V£yz).  We  note 
that  the  power  spectrums  for  eyy  or  eyz  are  obtained  from  that  of  exx  or  sxz  rotated  by  90°. 


Two  most  important  features  of  the  spectrums  are  their  frequency  bandwidth  and  rota¬ 
tional  symmetry.  If  the  sensor  responses  cannot  be  combined  to  form  a  rotationally  symmetric 
spectrum,  the  sensor  response  would  depend  on  stimulus  orientation  in  the  xy-plane,  which  is 
undesirable,  because  of  different  proportions  of  signal  to  noise,  (i.e.  assuming  that  the  noise  is 


53 


at  least  partially  independent  of  the  signal  amplitude). 


We  note  that  VElz  has  similar  shape  as  VPn  but  higher  bandwidth,  while  the  other  strains 
have  orientation  dependent  bandwidth.  This  excludes  candidate  (e)  since  for  VExy  (Fig.  3-5), 
frequencies  along  ±45°  are  attenuated.  Candidate  (f)  is  also  excluded  because  both  exx  and 
£yy  attenuate  frequencies  along  0  and  90°  and  using  both  will  therefore  not  produce  rotational 
symmetry.  It  is  worth  noting  that  the  directions  of  attenuation  coincide  with  the  conditions 
under  which  the  TFM  looses  rank  (see  section  3.2).  Candidate  (g)  is  therefore  the  only  one 
remaining  and  it  can  provide  rotational  symmetry  as  seen  in  Fig.  3-7.  The  importance  of 
having  full  rank  and  bandwidth  in  all  directions  becomes  apparent  when  one  considers  that 
loss  of  rank  means  that  infinitely  many  load  combinations  produce  identical  sensor  responses 
and  hence  make  inversion  impossible  for  those  cases,  and  lower  bandwidth  in  any  one  direction 
means  that  less  information  on  the  spatial  features  characterizing  that  dimension  of  the  object 
is  observable  by  the  sensors.  The  conclusion  is  therefore  that  given  the  assumption  that  we 
are  using  Cartesian  strain  sensors,  the  only  combination  that  works  for  incompressible  materials 
and  possesses  symmetry  is  ( pn ,  exz ,  eyz)  and  the  TFM  for  that  selection  is  (ref.  Appendix  A): 


Tsf  = 


(1  +  v)t 
EQ 


-zQ 


J^X 

{tt  -  iolz) 


JUy 


Q 

j  z 


L OxU!yZ  —  —  UJyz'j  jtlujyz 

1  T 

where  s  stands  for  a  normalized  sensor  signal  s  =  jgPn,  ?xz,  £yz\  ■ 


(3.4) 


Remark:  The  sensor  combination  (ozz,  ctx-,  ayz)  is  equivalent  to  ( pn ,  exz ,  eyz),  but  the  latter 
combination  can  be  considered  to  be  unique  and  optimal  because  in  almost  all  sensors,  it  is 
the  strains  that  are  directly  measured  and  stresses  are  then  inferred  from  the  strains.  The 
mean  normal  stress,  pn ,  is  however  needed  for  incompressible  materials  and  therefore  included. 
As  long  as  we  are  sensing  stresses  or  strains,  the  range  of  spatial  frequency  sensed  through 
each  component  is  of  the  same  order.  But  only  some  of  the  components  and  combinations  are 
rotationally  symmetric  about  the  2-axis. 


The  Structure  of  the  Sensor  TFM 


Analyzing  the  structure  of  the  TFM  in  Ecp  (3.4)  we  observe  that  the  e  nz  factor,  which  is 
common  to  all  the  terms,  expresses  how  the  signals  get  increasingly  lowpass  filtered  or  blurred 
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Pn 


with  depth.  It  is  a  consequence  of  the  model  including  the  resistance  of  the  medium  to  shear 
and  is  therefore  necessary  if  the  sensors  are  to  have  a  subsurface  location.  We  further  observe 
that  there  are  two  basic  variants  of  transfer  functions  in  the  sensor  transfer  function  matrix, 
which  we  will  call  T1  and  T2.  T1  is  characterized  by  its  norm  behaving  like  c~~lz ,  for  example 
TPiih  =  e~nz.  T2,  on  the  other  hand,  is  characterized  by  its  norm  behaving  like  |lU,|  . 

i  =  x ,  y,  e.g.  | TSxzfz\  =  \uix\  ze~Qz.  They  differ  in  that  the  T1  transfer  function  has  lowpass 
characteristics  but  the  T2  transfer  function  has  bandpass  characteristics.  In  the  spatial  domain 
this  means  that  the  response  of  the  T1  terms  is  spread  over  larger  area  since  a  predominantly 
low  frequency  content  indicates  slow  variations  in  the  spatial  domain.  The  response  of  the  T2 
terms  on  the  other  hand,  is  more  localized  and  they  respond  better  to  sudden  changes  in  the 
loading  since  they  have  higher  bandwidth. 
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Figure  3-3:  The  power  spectrum  of  the  x  normal  strain  V£xx.  This  power  spectrum  is  not 
rotationally  symmetric  as  it  behaves  like  VVn  in  the  ujx  =  ±uy  directions  but  as  Ve,,  in  the 
^.-direction  (along  ujy  =  0.) 


Terms  like  j—e 
Q 


-a. 


x,  y,  have  the  same  bandwidth  as  the  T1  and  but  terms  like 


n 


e  u-,  i,  k  =  x,  y  on  the  other  hand  are  of  type  T2. 


The  TFM  is  symmetric  with  respect  to  x  and  y,  since  the  halfspace  model  is  identical  for  x 
and  y.  Looking  at  the  TFM  linewise,  we  see  that  line  1,  ( pn  terms)  contains  only  T1  terms,  but 
lines  2  and  3  (sxz  and  eyz  terms)  contain  mostly  T2  terms.  T1  terms  appear  in  the  subdiagonal, 
in  T£izSt  and  T£yzfy.  These  T1  terms  serve  their  purpose  as  they  are  are  the  reason  that  the 
matrix  is  non-singular  for  0  =  0. 


Looking  at  the  TFM  columnwise,  and  column  3  first,  we  see  that  pn  term  in  line  1  is  type 
T1  and  the  sxz  and  eyz  terms  are  type  T2,  and  we  observe  that  they  are  proportional  to  the 
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Figure  3-4:  The  power  spectrum  of  the  2  normal  strain  Ve,z.  This  one  is  also  rotationally 
symmetric  as  VPn  (Fig  3-2)  but  it  has  higher  bandwidth. 

x  and  y  partial  derivatives  of  the  pn  term  (in  addition  also  multiplied  by  z).  The  same  is  true 
for  the  other  two  columns,  except  for  the  previously  mentioned  additional  T1  terms  on  the 
subdiagonal. 

In  summary,  the  TFM  can  be  interpreted  as  the  pn  terms  measuring  the  surface  force  as 
directly  as  possible  within  the  constraints  posed  by  the  subsurface  location  while  the  other  two 
sensor  signals  complement  the  pn  signal  by  measuring  its  x  and  y  spatial  derivatives. 
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Figure  3-5:  The  power  spectrum  of  the  xy  shear  strain  VZrxj.  It  is  not  rotationally  symmetric 
and  it  differs  from  the  other  spectrums  in  that  it  has  a  sharp  minimum  at  the  origin. 

3.3  The  Decoding  Solution 

The  Ideal  Case 

Ideally  the  solution  to  the  problem  would  be  obtained  by  simply  dividing  the  sensor  output  by 
the  transfer  function  in  the  frequency  domain,  i.e. 

00  (^£1  ^yi  fl)  —  T  z'js^CJxj 

where  0o  is  the  desired  surface  signal  that  needs  to  be  calculated  during  the  decoding  process 
(can  be  either  tractions,  displacements  or  a  combination),  T(u>x,u>y,  z)  is  the  transfer  function 
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Figure  3-6:  The  power  spectrum  of  the  xz  shear  strain  VExz.  This  power  spectrum  behaves 
like  V£z,  in  the  ^-direction,  but  as  VPn  in  the  o^-direction,  and  it  is  therefore  not  rotationally 
symmetric.  The  Vs  power  spectrum  can  be  obtained  from  this  one  by  a  permuationa  of  x,  y 
or  a  rotation  of  the  axis  by  ±90o- 

and  s  is  the  measured  signal.  As  can  be  seen  from  the  TFM  expressions,  the  transfer  function  is 
a  low-pass  filter  where  the  high-frequency  response  goes  exponentially  to  zero  as  Q  tends  to  oo. 
Therefore  the  inverted  transfer  function  will  become  arbitrarily  large  for  high  frequencies.  This 
is  very  impractical  and  will  lead  to  serious  errors  in  real  circumstances.  To  avoid  this  problem  we 
will  instead  use  regularized  solutions,  which  is  a  common  practice  in  treating  ill-posed  problems 
such  as  this  one.  Most  regularization  problems  are  formulated  as  scalar  problems,  here  however 
we  will  need  to  use  multivariate  regularization  which  happens  to  be  used  in  specialized  problems 
such  as  the  digital  restoration  of  multichannel  images  [21,  22], 
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Figure  3-7:  The  power  spectrum  of  the  sensor  signals  in  candidate  (g),  i.e.  'P(Pn<£xzt£  )  = 
T>Pn+'Pexz  +'P£y. .  We  observe  that  it  is  rotationally  symmetric  and  as  the  exz  and  eyz  components 
complement  each  other. 


The  Real  Case 


In  addition  to  the  features  of  the  ideal  problem,  the  real  problem  involves  the  unfavorable 
effects  of  noise  and  spatial  sampling.  Both  these  effects  can  be  effectively  minimized  using 
regularizing  agents  in  the  solution  as  shown  in  [80].  One  possible  way  of  solving  the  resulting 
integral  equations  is  to  discretize  the  integral  and  put  it  in  the  form  of  a  system  of  linear 
equations  [65,  63].  A  more  intuitive  way  is  to  construct  a  regularizing  operator  in  the  frequency 
domain  using  transfer  functions  as  follows.  The  solution  is  obtained  as  in  the  ideal  case  except 
for  the  multiplication  factor  ip(u>x,u)y,  z,a)  where  <x  is  a  regularizing  parameter,  i.e. 

^0  (ate  ,  CUy  )  'll)  (uq  ,  UJy  ,  ,  <7  )  T  [uJX1^ Uy  ,  z)  S  (tlTc  ,  CUy  ,  z)  . 
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where  is  an  estimate  of  the  surface  signal  vector  (inputs)  and  s  is  stimuli  (output)  vector. 

One  possible  family  of  regularizing  operators  that  can  be  constructed  are  weighted  norm 
filters  [80,  38]  which  are  of  the  form 


4>H,a(“x,iOy,z,a)  =  [T*T  +  aH}-1  T*T 

where  H(u)x,uy,  z)  is  weighting  function  to  be  determined  and  will  depend  on  the  problem  at 
hand,  the  sampling  interval  and  the  noise  characteristics.  The  solution  is  then 

<£o  =  [T*T  +  aH]~1T*s 

It  is  shown  in  [80]  and  in  further  detail  in  [38],  that  in  the  case  of  scalar  signals  this  regularizing 
agent  corresponds  to  the  minimization  of  a  functional.  This  generalizes  readily  in  the  vector 
case  to  the  functional 

M<t[oo,  s,  a]  =  [T(p o  —  s]r  [ T<f>o  —  s]  dcjxduy+ 

Cl  f_c o  f—QQ  ^y,  z)4>o(u>xi  LOy)dldxdu)y. 

This  means  that  in  the  solution  the  frequencies  of  the  solution  are  weighted  by  the  matrix  H. 


It  is  also  shown  in  [80,  38]  how  the  selection  of  the  parameter  a  and  the  matrix  H  can  be 
optimized  further.  If  we  assume  additive  sensor  noise  which  is  uncorrelated  with  the  solution, 
x,  y  and  load,  i.e.  s  =  sj  +  s jy  then  H  and  a  can  be  specified  such  the  expected  value  of  the 
squared  difference  is  minimized,  i.e. 


H  i|2 


where 


where  <?q  is  the  uncontaminated  signal.  The  solution  to  this  problem  is  the  multivariate  Wiener 
inverse  used  in  problems  which  appear  in  digital  restoration  of  multichannel  images.  It  is  shown 
in  [21,  22]  that  in  the  multivariate  case  it  is 

H(aJx,UJy,z)  —Ps  LOy  )Pn(wx,  LOy ,  z) 

where  Ps(ujx,u)y,  z)  is  the  power  spectrum  matrix  of  the  solution  and  Pn (ujx ,  uy )  is  the  power 
spectrum  matrix  of  the  sensor  noise.  The  multivariate  Wiener  inverse  then  finally  becomes 


C  =  [ PsT*T  +  Pn}-1  PsT*s. 


(3.5) 


The  need  for  the  power  spectrum  of  the  solution  indicates  that  the  solution  is  not  in  a  closed 
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form.  This  problem  can  be  circumvented  in  several  ways.  First,  the  power  spectrum  can  be  put 
equal  to  ql  where  I  is  the  identity  matrix  and  q  is  a  constant,  indicating  that  nothing  is  known 
about  the  solution  a  priori  and  the  solution  is  therefore  assumed  to  be  white  noise.  Second,  the 
solution  can  be  iterated.  Third,  all  prior  information  about  the  solution  that  can  be  included 
in  a  power  spectrum  can  be  used  there.  For  example,  suppose  Coulomb’s  law  of  friction  applies 
at  all  points  at  the  contact  interface  and  the  coefficient  of  friction  is  known.  The  solution  can 
then  be  constrained  to  satisfy  Coulomb’s  law  by  explicitly  expressing  the  correlation  between 
the  normal  and  the  tangential  tractions  in  the  components  of  the  Ps  matrix.  Let  us  further 
assume  that  it  is  known  that  the  object  is  sliding  over  the  surface  but  the  direction  of  motion 
in  the  (x,  y)-plane  is  unknown,  then  the  power  spectrum  matrix  would  be  of  the  form 


Ps  = 


p2f 

0 

pf 

0 

p) 

PS 

pf 

Pf 

1 

where  Pfz  is  the  power  spectrum  of  fz  and  pf  \s  the  friction  coefficient  between  the  object  and 
the  surface.  Using  this  Ps  in  Eq.  (3.5)  will  constrain  the  solution  to  satisfy  Coulomb’s  law,  i.e. 
fx  Pffz  and  fy  —  pjfz. 


Remark:  The  linear  filter  presented  here  is  based  on  the  assumption  that  the  noise  properties 
are  shift-invariant  and  uncorrelated  with  the  input  signal.  Also,  its  optimality  is  achieved  with 
respect  to  a  Euclidian  norm  functional.  The  shift-invariance  assumption  comes  naturally  with 
with  our  formulation  of  the  encoding  problem  as  a  contact  problem  on  the  homogeneous  semi- 
infinte  half-space.  However,  if  needed,  the  Wiener  filter  can  be  extended  to  correctly  treat 
cases  when  the  input  and  the  noise  are  correlated  [6],  Finally,  different  cost  functionals  can  be 
prescribed,  such  as  absolute  value,  p-norm  and  maximum  entropy  [10].  These  filters  should  be 
considered  when  the  Wiener  filter  has  proven  inadequate.  The  details  of  their  construction  is 
hence  considered  to  be  a  matter  of  implementation  and  beyond  the  scope  of  this  thesis. 


4.  Identification  of  Shape  and  the  Onset  of  Slip 


In  this  section  we  present  an  example  that  implements  the  theory  discussed  above  which  brings 
together  contact  mechanics  and  tactile  sensing.  The  simulations  have  two  main  objectives  i) 
shape  identifaction  and  ii)  the  detection  of  the  onset  of  slip.  In  both  cases  we  will  (a)  show  the 
necessity  and  usefulness  of  the  full  3D  formulation  presented  here,  especially  in  the  context  of 
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a  pseudo-dynamic  problem,  and  (b)  show  the  necessity  of  the  multivariate  regularizing  inverse 
in  the  solution  of  the  decoding  problem. 

The  example  problem  chosen  is  that  of  the  incipient  relative  sliding  of  two  elastic  bodies  in 
contact.  This  problem  is  particularily  interesting  in  the  context  of  haptics,  as  the  identification 
of  object  shape  and  the  prediction  of  the  onset  of  and  the  prevention  of  slip  is  a  fundamental 
problem  for  both  human  and  robotic  exploration  and  manipulation.  It  has  also  been  observed 
by  [66,  67]  that  objects  of  general  shape  will  have  qualitatively  similar  behavior  as  spherical 
objects,  so  the  results  presented  here  will  extend  to  general  objects. 


The  problem  demonstrated  here  is  the  indentation  of  an  infinite  half  space  by  an  object 
shaped  as  a  paraboloid,  with  simultaneous  tangential  (Fxo,Fyo)  and  normal  (FZQ)  load.  Assum¬ 
ing  Coulomb  friction  with  friction  coefficient  /ij.  the  condition  for  local  slip  can  be  expressed 
as 

f  <  y/fx(xFj)  +  fy(x,y)  =»  no  slip 
Vffz(x,y)l  _  (3.6) 

l  >  \Jfl{x,y)  +  fy(x,y)  =>  Slip  occurs 

It  is  shown  in  [36]  that  the  unique  solution  to  this  problem  is  that  local  slip  occurs  on  the  outer 
edge  of  the  contact  area  where  c  <  r  <  a,  (r  =  \/x2  y 2)  where 


c 


a 


(3.7) 


is  the  relative  size  of  the  non-slip  region  and  no  slip  occurs  for  r  <  c,  (see  Fig  3-13).  The 
resulting  surface  tractions  are 


.fz{x,  y) 


ft(x,y) 


for  0  <  ?•  <  a 


where  c  <  r  <  a 


where 


r  <  c 


where  ft  =  f.  —  0  for  r  >  a,  and  ft(x,y)  is  directed  along  [Fxo ,  Fyo]r  with  ft(x,y )  = 
\ffxix,y)  +  fy(x,y)- 
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4.1  Identification  of  shape 


We  use  here  the  optimal  sensor  combination  obtained  in  section  3.2  and  hence  the  encoding 
part  consists  of  predicting  the  sensor  signal  set  (^pn,  syz)  by  using  the  traction  vector 
/  above  with  the  TFM,  Tsf  given  in  Eq.  (3.4).  The  form  of  these  sensor  signals  is  shown  in 
Fig.  3-8(a)-(c).  The  decoding  part  consists  of  reconstructing  the  traction  and  the  displacement 
vectors,  from  noise  contaminated  sensor  signals.  The  displacements  will  give  us  the  shape  of 
the  object  where  it  is  in  contact  with  the  surface,  but  that  is  by  definition  where  the  traction  is 
non-zero,  and  hence  both  displacements  and  tractions  are  needed  for  object  shape  identification. 

We  will  now  illustrate  the  necessity  of  regularization  in  the  decoding  algorithm  by  adding 
noise  to  the  sensor  signals.  In  this  example  we  use  a  high-pass  filtered  Gaussian  noise,  and 
in  Figs.  3-8(d)-(e)  we  show  the  resulting  noise  contaminated  sensor  signals.  Profiles  of  an 
unregularized  solution  obtained  by  inverting  Tsf  in  Eq.  (3.4)  is  shown  in  Fig.  3-10  and  profiles 
of  a  regularized  solution  obtained  using  the  Wiener  inverse  given  in  Eq.  (3.5)  is  shown  in  Fig. 
3-9.  Surface  plots  of  the  same  solutions  are  shown  in  Figs.  3-11  and  3-12.  The  displacements 
are  shown  in  Fig.  3-11  and  the  surface  tractions  in  Fig.  3-12.  We  observe  that  the  regu¬ 
larized  inverse  is  quite  close  to  the  uncontaminated  solution  but  the  unregularized  one  is  not 
recognizable. 


4.2  Onset  of  Slip 


As  the  total  tangential  load,  Ft0  =  F'“0  +  F~0,  takes  increasing  values,  the  size  of  the  non-slip 

region  decreases  and  when  Fto  =  pjFZo  global  slip  will  occur.  Therefore,  knowing  the  relative 
size  of  the  non-slip  region  will  enable  us  to  estimate  how  close  a  particular  loading  situation 
is  to  slipping.  The  size  of  the  non-slip  region  can  be  inferred  from  the  tangential  traction 
distribution  as  there  will  a  sharp  edge  along  r  =  c  in  that  distribution.  This  edge  can  be  seen 
as  two  peaks  on  the  2D  profile  plot  of  ft{x,y)  in  Fig.  3-9  and  as  an  edge  in  the  /f-surface  plot 
in  Fig.  3-12.  We  can  measure  a  and  c  from  the  reconstructed  /t(x,  y )  with  the  help  of  standard 
edge  detecting  techniques.  That  information  along  with  the  total  tangential  and  normal  loads, 
Fto ,  F-0,  can  either  be  used  in  Eq.  3.7  to  solve  for  /z/  or  compared  to  the  line  plot  in  Fig.  3-13 
to  estimate  how  close  the  contact  is  to  global  slip.  For  example  we  see  from  Fig.  3-13,  that  if 
c  =  0.5a  then  the  tangential  load  has  reached  87.5%  of  its  maximum  value  (Fto  =  0.875 y.fF:o). 
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(a)  Noise  free  sensor  signal  pn  (b)  Noise  free  sensor  signal  £xz  (c)  Noise  free  sensor  signal  evz 


Figure  3-8:  Sensors  signals  without  (a)  -  (c)  and  with  (d)  -  (f)  noise. 

5.  Conclusions 


In  this  chapter  we  have  delineated  a  basic  computational  theory  of  haptics  that  is  common 
to  humans  and  robots.  The  tactile  identification  problem  is  separated  into  the  encoding  and 
decoding  problems  and  solutions  for  both  are  given.  Analysis  of  the  results  is  shown  to  lead 
to  a  unique  and  optimal  combination  of  sensors  suitable  for  tactile  sensing.  An  optimized 
multivariate  regularizing  algorithm  for  the  solution  of  the  decoding  problem  is  also  presented. 
Using  simulations,  it  is  shown  how  the  formulation  used  here  can  be  applied  to  predict  the 
onset  of  slip  using  tactile  sensors. 

The  conclusions  of  this  chapter  can  thus  be  summarized  as  follows: 


1.  To  sucessfully  decode  a  general  tactile  3D  signal  from  a  tactile  sensor  array,  a  mini¬ 
mum  of  three  independent  sensors  are  needed;  ( pn ,  £xz,  eyz)  is  a  unique  combination  of 
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fx  along  y=0 


x 


ux  along  y-0 


x 


fz  along  y=0 


uz  along  y=0 


Figure  3-9:  Shows  profiles  of  solutions  /x,  fz ,  ux,  uz  along  y  =  0  to  the  decoding  problem 
obtained  using  an  unregularized  inverse.  We  observe  that  the  tractions  are  hardly  recognizable 
although  the  solutions  for  the  displacements,  ux,  uz  are  closer  to  the  true  solution  fx,  fz. 


stress/strain  sensory  signals  which  do  not  have  a  directional  bias. 

2.  The  transfer  function  approach  developed  here  allows  us  to  look  at  the  qualitative  prop¬ 
erties  and  behavior  of  the  solution,  irrespective  of  the  exact  form  of  the  load  or  object 
shape. 

3.  The  two  items  above  lead  to  the  transfer  function  matrix  (TFM)  formulation,  a  convenient 
and  compact  formulation  that  is  well  known  in  control  theory.  The  decoding  problem  is 
then  reduced,  essentially,  to  the  inversion  of  the  (TFM). 

4.  Since  the  inversion  is  ill-posed  in  the  presence  of  noise  and  discrete  sampling,  regulariza¬ 
tion  is  needed,  for  which  the  multivariate  Wiener  inverse  exists  as  a  standard  tool. 

5.  Simulations  presented  in  section  4  show  that  the  above  can  be  used  to  calculate  surface 
shape,  load  distribution,  contact  area,  and  further  has  a  possible  application  in  estimating 
the  onset  of  slip. 
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f x  along  y=0 


x 

ux  along  y=0 


x 


fz  along  y=0 


ut  along  y=0 


Figure  3-10:  Shows  profiles  of  solutions  to  the  decoding  problem  obtained  using  the  regularized 
Wiener  inverse.  We  observe  that  it  coincides  almost  everywhere  with  the  true  solution  and  that 
the  error  is  therefore  small. 
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(a)  True  ux  displacement 


(b)  Estimated  ux  (Wiener)  (c)  Estimated  ux  (unregularized) 


u  0 
>'  -5  -5  x 


Figure  3-11:  Shows  surface  plots  of  decoded  surface  displacements  ux  in  (a)-(c),  and  uz  in  (d)- 
(f).  As  in  Figs.  3-10  and  3-9  we  observe  that  the  true  solution  and  the  solution  obtained  using 
the  Wiener  inverse  look  very  similar  but  the  solution  obtained  using  an  unregularized  inverse 
is  considerably  worse. 
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(a)  True  fx  loading  (b)  Estimated  fx  (Wiener)  (c)  Estimated/^  (unregularized) 


y  -5  -5  x  y  -5  -5  X  >’  -5  -5  x 

Figure  3-12:  Shows  surface  plots  of  decoded  surface  load  distributions  fx  in  (a)-(c),  and  fz  in 
(d)-(f).  Again  observe  that  the  true  solution  and  the  solution  obtained  using  the  Wiener  inverse 
look  very  similar  but  the  solution  obtained  using  an  unregularized  inverse  is  mostly  noise. 
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Relative  Slip  Area  vs.  Relative  Tangential  Load 


Figure  3-13:  The  relation  between  the  relative  tangential  load  FtJ(nfFZ0),  where  Fto  =  ( F £  + 
Fy  )1//2,  and  the  relative  no-slip  radius  c/a  =  1  —  F/0/(jj.fFZf))'i.  As  an  example,  this  relation 
indicates  that  when  the  no-slip  radius  is  equal  to  half  the  contact  radius  (c/a  =  1/2)  then 
Ft0  =  0.875 H;FZ0  or  that  the  tangential  load  is  87.5%  of  the  maximum  tangential  load. 
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Chapter  4 


Dynamic  Identification 


1.  Introduction 


An  understanding  on  the  nature  of  touch,  especially  the  mechanics  and  the  information  transfer 
involved,  is  important  for  both  physiological  and  robotics  research.  The  workings  of  human 
touch  are  presently  not  fully  known  and  a  computational  model  providing  theoretical  basis 
would  be  useful  when  interpreting  experimental  results  as  a  tool  to  generate  and  validate  con¬ 
sistent  hypothesis.  In  robotics,  a  theoretical  model  would  be  of  use  in  the  design  and  devel¬ 
opment  of  sensor  arrays,  addition  to  providing  data  processing  algorithms  and  paradigms  for 
their  implementation.  Finally,  in  the  emerging  fields  of  virtual  environments  and  telerobotics 
depend  considerably  on  detailed  knowledge  on  the  nature  of  the  interaction  of  humans  and 
machines  through  touch.  The  above  leads  us  to  conclude  that  a  computational  theory  which 
would  cover  the  common  aspects  of  humans  and  robotic  touch  could  be  of  considerable  and  far 
reaching  value  and  we  will  in  this  chapter  discuss  the  details  of  such  a  theory  which  relates  to 
tactile  sensing  under  dynamic  contact  conditions. 

In  chapter  3,  we  have  treated  a  related  problem,  where  we  assumed  that  the  loading  is  static 
or  quasi-static.  The  dynamic  effects  become  important  when  the  loading  varies  significantly  on 
a  timescale  measured  in  units  of  time  constants  characteristic  of  the  medium  in  the  fingerpad. 
In  this  chapter  we  will  include  the  dynamic  effects  of  inertia  and  viscoelasticity  in  our  analysis 
of  the  tactile  sensing  problem.  We  will  then  delineate  the  cases  in,  and  conditions  under  which 


these  effects  become  significant  and  when  the  can  be  neglected  -  giving  examples  where  the 
static  solution  given  in  chapter  3  is  no  longer  valid  and  this  dynamic  formulation  becomes 
necessary. 


We  use  the  same  framework  and  definition  of  terms  as  in  chapter  3,  which  we  state  here 
again  for  review:  Two  stages  in  tactual  information  processing,  encoding  and  decoding  make 
up  the  identification  problem ,  which  is  so  called  because  it  involves  the  determination  of  the 
state  of  the  object  and  its  relationship  to  the  tactual  sensory  system.  The  system  discussed 
here  differs  in  some  ways  from  traditional  dynamic  systems.  The  most  prominent  difference  is 
that  it  has  both  spatial  and  temporal  variables,  whereas  traditional  lumped  parameter  systems 
only  have  temporal  ones.  In  both  natural  and  man-made  tactile  systems  the  encoding  process 
is  the  transduction  of  the  mechanical  stimulus  into  electric  signals  that  contain,  in  a  coded 
form,  information  about  the  stimulus  features.  The  inverse  problem,  i.e.  the  decoding  process 
or  the  calculation  of  the  surface  shape  or  surface  tractions  from  the  transduced  signals  is  in 
many  respects  significantly  more  difficult.  The  encoding  is  essentially  a  spatial  low-pass  filtering 
process  which  indicates  that  the  decoding  will  be  ill-posed  in  the  sense  that  the  higher  frequency 
components  of  the  signals  will  be  prone  to  noise  contamination  and  special  methods  have  to  be 
applied  for  successful  decoding.  This  applies  equally  to  the  static  formulation  given  in  chapter 
3,  but  several  features  are  added  to  this  problem  with  the  dynamic  formulation.  Temporal 
properties  such  as  inertia  and  viscosity  now  play  a  role  in  the  dynamic  encoding  solution,  and 
the  decoding  problem  also  involves  the  causality  constraint  on  the  temporal  dimension  that 
needs  to  treated  accordingly. 

Several  papers  can  be  found  in  the  robotics  literature  on  the  mechanistic  modeling,  encoding 
and  decoding  of  surface  conditions  from  subsurface  measurements  using  tactile  sensors  under 
static  conditions  and  the  reader  is  refered  to  chapter  3  for  an  overview.  Papers  where  dynamic 
effects  are  included  are  much  more  scarce,  but  in  [28]  the  encoding  and  decoding  of  a  stress  rate 
sensor  are  investigated  and  a  scalar  Wiener  inverse  is  suggested  for  the  decoding.  In  [84]  the 
response  of  a  2D  plane  strain  model  to  vibrating  sinusoidal  gratings  is  analyzed  and  it  compared 
to  the  sensitivity  of  the  human  fingerpad  at  different  vibration  and  grating  frequencies. 

All  previously  published  results  on  the  reconstruction  of  general  time- varying  surface  loads 
and  profiles  from  subsurface  tactile  signals  do  not  solve  the  full  3D  problem  including  the  effects 
of  inertia  and  viscoelasticity  and  the  effects  and  interaction  of  shear  loads.  In  this  chapter  we 
will  address  those  aspects  of  the  mechanistic  modeling,  encoding,  and  the  decoding  problems, 
and  the  contributions  of  this  chapter  can  be  itemized  as  follows: 
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1.  A  complete  solution  to  the  3D  half-space  problem  is  given  for  arbitrary  3D  dynamic 
loading  conditions  and  any  linearly  elastic  material.  The  solution  is  formulated  using 
transfer  function  matrices  in  the  domains  of  spatial  and  temporal  frequencies  which  allows 
us  to  look  at  the  qualitative  properties  and  behavior  of  the  solution,  irrespective  of  the 
exact  form  of  the  load  or  object  shape. 

2.  The  effects  of  inertia  and  viscoelasticity  are  analyzed  and  it  is  shown  that  inertial  effects 
can  be  neglected  whereas  the  visocelastic  effects  are  dominant  for  low  temporal  frequen¬ 
cies.  The  use  of  linear  viscoelastic  models  further  leads  to  solutions  that  can  be  expressed 
using  Transfer  Function  Matrices  that  are  rational  in  the  temporal  frequency. 

3.  For  decoding,  we  can  then  apply  the  extensively  developed  and  numerically  efficient 
Kalman  filter  to  estimate  the  inputs  which  at  the  same  time  avoids  the  problems  as¬ 
sociated  with  the  ill-posedness  caused  by  the  spatial  low-pass  filtering. 


In  Section  2,  we  formulate  and  solve  the  3D  half-space  elastic  model,  which  serves  as  a 
description  of  the  encoding  process.  For  the  sake  of  clarity,  the  mathematical  expressions  are 
detailed  in  Appendix  C.  We  analyze  the  solutions  in  order  to  assess  if  and  when  the  effects 
of  inertia  can  be  neglected  and  further  show  how  viscoelastic  effects  are  included.  In  Section 
3,  we  solve  and  analyze  the  decoding  problem,  and  in  Section  4,  we  give  example  numerical 
simulations  of  the  methods  proposed  in  Section  3.  Conclusions  are  given  in  section  5. 


2.  The  Encoding  Problem 


From  theoretical  viewpoint,  tactile  sensor  response  to  any  given  stimulus  can  be  predicted  if 
we  have  the  following: 

1.  A  mechanistic  model  of  the  deformable  medium  which  can  be  used  to  calculate  reliably 
the  stresses  and  strains  at  each  instant  and  at  each  point  in  the  medium  for  a  given  time 
history  of  the  mechanical  stimulus. 

2.  A  model  of  the  sensor  that  provides  the  relationship  between  the  relevant  stimulus,  i.e. ,  a 
particular  combination  of  stresses  and  strains  in  the  local  neighborhood  of  a  sensor  that 
it  is  responsive  to  and  the  temporal  characteristics  of  the  sensor  response. 
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Figure  4-1:  Shows  the  semi-infinite  halfspace  with  surface  tractions  and  subsurface  sensors. 
The  x-  and  //-dimensions  extend  to  infinity  in  both  negative  and  positive  directions.  The  z- 
dimension  extends  from  0  to  positive  infinity.  The  tractions  are  distributed  over  the  contact 
area,  each  point  being  subjected  to  3-dimensional  loading  f(x,y,t)  =  [fx.  fy ,  fz]T . 


For  example,  an  appropriate  mechanistic  model  of  the  human  skin  and  the  soft  tissues  under¬ 
neath  enables  us  to  investigate  one  of  the  important  questions  in  cutaneous  neurophysiology, 
that  of  identifying  the  relevant  stimulus  that  causes  each  receptor  type  to  respond  [35,  75], 
where  only  the  temporal  characteristics  are  known  but  the  particular  combination  of  strains  is 
not  agreed  upon. 


In  this  section  we  will  find  and  analyze  the  dynamic  response  of  the  elastic  halfspace  sub¬ 
jected  to  a  time-varying  load  distribution,  and  based  on  that  give  conditions  under  which  we 
can  neglect  the  effects  of  inertia  and  how  viscoelasticity  can  be  included  in  the  solution. 


The  problem  is  set  up  so  that  the  boundary  surface  coincides  with  the  xy-plane  and  the 
positive  r-axis  points  into  the  medium  (See  Fig.  4-1).  The  material  is  assumed  to  be  linearly 
elastic  and  homogeneous.  It  is  also  assumed  to  extend  to  infinity  in  both  negative  and  positive 
x.  y  directions  and  in  the  positive  2  direction.  The  surface  load  is  described  by  the  time- varying 
traction  vector  field  f(x,y,t.)  —  [ fx{x,y,t ),  fy(x,y,t ),  fz(x,  y,  t)]T  and  the  surface  displacement 
field  as  uo(x,  y ,  t)  =  [ux(.r,  y,  0,  t),  uy(x ,  y,  0,  t),  uz(x,  y,  0,  t)]T .  In  our  formulation,  the  boundary 
conditions  can  be  forces,  displacements  or  mixed,  i.e.  when  any  3  elements  from  both  /  and  «o 
are  given.  As  in  chapter  3  we  intend  to  present  the  solution  in  the  space  of  spatial  frequencies 
and  use  transfer  functions  to  simplify  solutions  and  gain  further  insight. 


2.1  Transfer  Function  Matrices  (TFM) 


We  solve  the  encoding  problem  by  transforming  the  fundamental  differential  equations  of  equi¬ 
librium  and  the  stress-strain  relationship  equations  (Generalized  Hooke’s  law)  using  the  Fourier- 
Laplace  transform  over  ( x,  y,t )  space.  This  is  analogous  to  what  was  done  to  solve  the  static 
problem  in  chapter  3  and  examples  of  this  approach  can  also  be  seen  in  [16]  and  [71]. 

The  differential  equations  of  equilibrium1  are  given  by 

Oij,j  =  pui  i,j  —  x,  y,  z 

(4.1) 

and  the  linear  stress  strain  relations  according  to  Hooke’s  law  are 

^ ij  —  Ac ij 6 ij  +  2 fi £ jj  i,j  —  x 

(4.2) 

where 

1  ,  ... 

Sij  —  ^  (tijj  +  Ujj)  —  x ,  y,  z 

(4.3) 

and  Sij  —  1  if  and  only  if  i  —  j,  and  is  zero  otherwise.  Further,  rho  is  density  and  the  Lame 
constants  A  and  y  (shear  modulus)  are  in  terms  of  Young’s  modulus  and  Poisson  ratio: 


x  vE  E 

X~  (1  +  n)(l  —  2v) '  2(1  +  1/)' 

The  details  of  the  derivation  are  given  in  Appendix  C  for  the  purely  elastic  case.  The 
complete  solution  when  viscoelastic  effects  are  also  present  is  treated  in  section  2.3. 

As  is  done  in  chapter  3,  we  employ  here  the  concept  of  a  transfer  function  matrix  (TFM).  The 
transfer  function  matrices  we  give  here  have  the  Fourier-Laplace  transform  of  the  components 
of  either  f(x,y,t )  or  uo(x,y,t)  as  the  inputs  and  any  displacements,  strains  or  stresses  in  the 
medium  as  the  outputs.  Each  entry  of  the  TFM  would  then  correspond  to  the  contribution  of 
a  particular  input  component  to  the  corresponding  output  component.  For  example  the  entry 
in  the  TFM  for  the  ^-displacement  resulting  from  the  r-load,  or  u*(x,y,  z,t)  we  would  write 

^yi  T  ^')  =  (tnx,  <Xy ,  z ,  s)/a.(tnX)  k-’y,  s) 

1Here  we  use  tensor  notation  for  compactness,  where  a  ”,  j"  in  subscript  means  a  partial  derivative  with 
respect  to  the  independent  variable  represented  by  j  and  repeated  symbols  in  subscript  stand  for  summation 
over  all  indices. 
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to  obtain  the  Fourier-Laplace  transform  uxz  (u>x,  uiy,  z,  s)  of  ux(x,y,z,t)  from  the  transfer  func¬ 
tion  Tuxf(ux,uy,z,s)  and  the  Fourier-Laplace  transform  fx(u:x,u)y,s)  of  the  surface  shear  load 
fx(x,  y,  t).  Here,  ux  and  uiy  correspond  to  spatial  frequencies  in  the  x  and  y  directions  respec¬ 
tively,  and  s  is  the  Laplace  variable  (related  to  the  temporal  frequency  by  s  —  ju>).  Then  we 
can  write 


u  =  Tuf  f  £  =  Teff  a  =  Tafj  pn  =  TPnJf 

u  =  TUUouo  £  =  T£Uquo  (j  =  Tauov,Q  pn  =  Tpn1lou  o 


(4.4) 


where  the  subscript  /  indicates  that  surface  tractions  are  inputs  and  the  subscript  uo  indicates 
that  surface  displacements  are  inputs. 


The  equations  in  (4.4)  represent  the  solution  for  the  encoding  problem,  because  given 
f(x,y,t),  uo(x,y,t)  or  a  combination  of  3  elements  from  the  two,  we  can  pi'edict  u(x,y,z,t), 
e(x,y,z,t),  a(x,y,z,t )  and  pn(x,y,  z,t).  If  the  relevant  stimulus,  i.e.  a  particular  combination 
of  stresses  and  strains  in  the  local  neighborhood  of  the  sensor  that  it  is  responsive  to,  is  known, 
we  can  predict  the  sensor  response.  For  example,  if  the  sensor  response  is  linearly  related  to  the 
strain  ezz,  then  it  is  the  relevant  stimulus  for  that  sensor.  On  the  other  hand,  if  the  relevant 
stimulus  needs  to  be  identified,  we  can  use  Eq.  (4.4)  to  generate  hypotheses  that  can  be  tested 
against  the  empirically  determined  sensor  response,  as  is  done  in  [35,  75], 


2.2  Analysis  of  the  Encoding  Solution 


The  above  solution  to  the  decoding  problem  gives  us  every  component  of  both  the  strain  and 
stress  tensors  but  since  the  input  is  a  3-dimensional  vector  we  only  need  to  select  three  signals, 
or  three  independent  combinations  of  signal,  which  the  sensors  should  transduce.  In  chapter  3 
the  selection  of  suitable  sensors  was  discussed  and  there  we  showed  using  symmetry  arguments 
that  the  sensors  should  optimally  transduce  three  signals:  the  mean  normal  pressure,  pn  and  two 
~  shear  strains  sxz  and  eyz.  The  arguments  made  in  chapter  3  will  also  hold  for  the  dynamic 
encoding  solution  as  introduction  of  time  as  a  variable  does  not  destroy  spatial  symmetry2. 
We  therefore  use  transfer  functions  for  the  signals  above  as  examples  in  our  analysis  of  the 
encoding  solution.  We  however  note  that  the  same  analysis  could  be  applied  to  all  the  other 
strain  and  stress  signals  and  the  conclusions  would  be  the  same.  Hence,  we  will  analyze  two 
transfer  functions,  one  between  the  mean  normal  stress  pn  and  the  normal  load  /.,  and  the 

'The  Laplace  variable  s  only  appears  in  the  sensor  signals  as  some  constant  times  s 2  being  added  to  f 22  and 
will  therefore  not  introduce  any  spatial  asymmetry 
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other  between  the  22-shear  strain,  exz  and  fz. 


Transfer  Functions:  In  the  dynamic  solution  the  transfer  functions  for  those  two  variables 
become 


LPnfz  ~ 


£xz  f z 


3 02  —  4  s2  (fi2  +  5Csh2s2)e  ”12 

1  joJxni( ft2  +  jcJ^s2)  (e~n2Z  -  e~niz) 
F3 


(4.5) 


where 


F3  =  (fi2  +  \cjs2)2  -  f l2nxn2, 


n  1 


=  \j$l2  4-  cp2s2, 


n2 


\/tt2+cJs2, 


and  cp  =  t/(A  +  2p)/ p  and  cs^  =  \/p/ p  are  the  dilational-  and  shear-wave  wavespeeds  (see 
[81]),  with  A  and  fi  as  bulk  and  shear  modulus  respectively,  p  is  the  density,  and  0  =  cp/csh- 
As  before  ux  and  uoy  are  the  spatial  frequencies,  <4  =  +  tu2,  and  s  is  the  Laplace  transform 

variable,  related  to  the  temporal  frequency  by  s  =  jw.  The  solution  obtained  here  should 
therefore  approach  the  static  solution  when  s  — »  0,  and  we  show  that  does  so  in  Appendix  D. 


Analysis:  For  analysis  of  the  transfer  function,  it  is  helpful  to  normalize  and  measure  distances 
and  displacements  in  units  of  shear-wave  wavelengths.  Let  s  =  ju,  A sh  =  csh/u>  and  hence 


Z  —  zX  sh  —  ZCshjuJ 

—  Cl  j A  sh  —  Clid  /  Cgfi 


M x  —  &x/^sh  —  &x^/c$h 

The  temporal  frequency,  w,  can  now  be  canceled  out  from  the  transfer  functions,  which  then 
become 


lPnfz 


3/32  —  4  (n2-|)e-"12 

602  [(^2  _  1/2)2  _  Q2^lA2j 


l  jLJxfll(Cl2  -  i)  | 

[e-*2*  -  e“ 

2/r 

\Cl2  -  1/2)' 

2  —  Cl2fi\n2 

with  h\  =  (Cl2  —  0  2 )1'/2,  and  ?i2  =  (Cl2  —  1)1|/2. 


We  observe  that  all  the  terms  in  the  transfer  function  are  real  for  f2  >  1.  This  can  be 
interpreted  as  those  spatial  frequencies  which  satisfy  Clcsh  >  u>  that  are  not  involved  in  the 
transfer  of  energy  away  from  the  source  and  that  the  halfspace  behaves  as  an  elastic  spring  for 
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Figure  4-2:  Shows  the  modeshape  for  the  lowest  eigenfrequencies  of  a  domain  with  length  L,  assuming 
either  fixed  or  free  boundary  conditions  at  x  =  ±L/2.  The  maximum  wavelength  is  Amax  =  2 L  and 
hence  the  lowest  non-zero  spatial  eigenfrequency  is  fim;n  =  tr/L. 


those  frequencies. 


The  denominator  of  the  transfer  functions,  however  vanishes  for  Cl  =  fi r,  something  we 
know  because  the  denominator  is  the  characteristic  equation  for  the  Rayleigh  wavespeed  (see 
[81]).  The  transfer  functions  therefore  become  singular  and  a  resonance  phenomena  occurs 
when  Cl  =  cr/csh,  where  cr  is  the  Rayleigh3  wavespeed. 


Inertia  and  Finite  Domains:  The  above  can  be  used  to  establish  conditions  under  which 
one  can  neglect  inertia,  especially  in  cases  where  the  solution  in  Eq.  (4.5)  is  applied  to  finite 
domains.  As  an  example,  for  a  rectangular  domain,  i.e.  (x,y)  €  [—L/2,  L/2]  x  [— L/2,  L/2], 
only  a  discrete  set  of  spatial  frequencies  will  give  valid  solutions.  The  lowest  nonzero  spatial 
frequency  is  of  particular  interest  to  us,  and  we  have  plotted  corresponding  the  modeshapes 
assuming  either  free  or  fixed  boundary  conditions  at  x,  y  =  ±L/ 2  in  Fig.  4-2. 

We  see  that  the  longest  possible4  wavelength  is  Amax  =  2 L  which  is  equivalent  to  f2m;n  = 

3The  value  of  cr  depends  on  3  =  cp/csh .  For  materials  that  are  nearly  incompressible  (f3  — >  co)  we  have  that 
cr  «  0.95csh. 

4This  is  assuming  fixed  or  free  boundary  conditions.  A  third  kind  of  boundary  conditions,  elastic  boundary 
condition  is  also  possible,  and  the  lowest  eigenfrequency  will  then  depend  on  the  stiffness  of  the  boundary 


7 t/L.  We  therefore  have  that  f2  <  cn/csh  is  satisfied  if 

.  VCR 

LO  <  Ucr  =  — (4.6) 
giving  us  a  condition  as  to  when  we  can  ignore  the  effects  of  inertia  for  the  semi-infinite  halfspace. 

Figures  4-3  and  4-4  show  how  the  condition  above  appears  in  a  Bode  plot  of  the  transfer 
functions  in  Ecp  (4.5).  The  numerical  values  were  chosen  so  that  the  system  resembled  the 
human  fingerpad  in  size  and  material  properties,  i.e.  L  =  12.8  [mm],  p  =  2.8- 104  Pa,  v  —  4.995, 
and  density  p  —  1000  [kg/m3],  which  substituted  into  the  condition  in  Ecp  (4.6)  gives  that 
luct  —  1095  rads-1  or  that  inertia  can  be  ignored  and  the  quasi-static  solution  used  if  the 
bandwidth  of  the  loading  is  less  than  174  Hz.  We  observe  in  Figs.  4-3  and  4-4  that  the  static 
and  the  dynamic  solution  are  identical  below  the  temporal  frequency  predicted,  but  that  there 
are  complex  pole-zero  pairs  (i.e.  (s2  +  z2)/(s 2  +  p2))  which  are  very  close  to  each  other  at  the 
critical  frequency.  Viscoelastic  behavior  for  frequencies  lower  than  that  critical  frequency  can 
now  be  modeled  by  letting  the  elastic  constants  in  the  static  solution,  given  in  chapter  3,  be  rate 
dependent,  which  will  give  some  shape  to  the  transfer  functions  in  that  range,  characterizing  the 
particular  material  being  modeled.  We  will  see  in  section  3  that  this  will  simplify  the  decoding 
part  significantly  and  in  fact  make  it  numerically  feasible. 


2.3  Viscoelasticity 


Theory  of  viscoelasticity  extends  the  theory  of  elasticity  to  include  cases  where  the  stress-strain 
relation  behavior  is  not  time  independent.  It  all  boils  down  to  the  choice  of  timescale,  and 
for  shorter  intervals  of  interest,  it  may  be  necessary  to  include  the  rate  dependent  behavior 
of  materials.  The  description  linear  viscoelastic  material  applies  to  materials  for  which  the 
rate  dependent  behavior  can  be  adequately  described  explained  by  expressing  the  stress-strain 
relation  in  the  viscoelastic  form  of  Hooke’s  law, 

a(s)  =  Q{s)e(s). 

Where  a(s )  and  e(s)  are  the  Laplace  transforms  of  the  stress  and  the  strain,  respectively.  Q(s) 
is  called  the  Operational  relaxance  and  alternatively  Q(s)/s  is  called  the  Operational  impedance. 
The  main  advantage  in  using  the  viscoelastic  form  of  Hooke’s  law  is  that  we  then  can  apply  the 


condition.  There  will  be  a  minimum  eigenvalue  similar  to  the  other  boundary  conditions,  with  the  exception  of 
the  singular  case  when  the  stiffness  of  the  boundary  condition  matches  that  of  the  material. 


Mean  Normal  Pressure  pn  static  —  dynamic 


T 


100  h 


o 

“O 


-100 


10' 


10 


/  [tfc] 


10 


x 


Figure  4-3:  Shows  a  Bode  plot  of  the  mean  normal  pressure,  pn  for  both  the  static  and  dynamic  solution 
at  z  =  1  mm  below  the  loading  point.  The  amplitude  is  normalized  with  the  static  response.  We  observe 
that  the  solutions  are  identical  up  to  a  critical  frequency,  given  by  u)cr  =  tt cr/L.  (  /  <  174  Hz  when 
estimates  of  the  material  parameters  for  the  human  fingerpad  are  used. 


correspondence  principle.  According  to  the  principle,  if  an  elastic  solution  to  a  stress  analysis 
problem  is  known,  substitution  of  the  appropriate  Laplace  transforms  for  the  material  constants 
employed  in  the  elastic  solution  furnishes  the  viscoelastic  solutions  [83]. 


Hence,  the  effects  of  viscoelasticity  can  be  incorporated  into  the  elastodynamic  solution  in 
Appendix  C,  simply  by  allowing  the  material  constants  to  be  frequency  dependent,  i.e.  p  and  0 
become  p(s)  and  0(s).  Similarly  the  effects  of  viscoelasticity  can  included  in  the  static  solution 
given  in  chapter  3. 


The  material  constants  are  in  general  both  independently  rate  independent,  p  describing 
shear  motion  and  0  the  dilational  behavior.  It  is  however  often  assumed  that  only  p  is  rate 
dependent.  This  assumption  is  reasonable  for  materials  that  are  near  being  incompressible, 
(v  >  0.4)  which  is  the  case  for  polymers  and  other  materials  we  intent  to  apply  this  theory  to. 
The  simplest  model  that  captures  the  behavior  of  a  linearly  viscoelastic  solid  is  the  3-parameter 


80 


5 


xz- Shear  Strain  exz  static  —  dynamic 

1  1  1  I - - - - - * - 1 — '  '  '<  '  l - - - - 


-L 


/  [/fc] 


10 


J. 


Figure  4-4:  Shows  a  Bode  plot  for  the  :rz-shear  strain,  exz  for  both  the  dynamic  and  static  solutions. 
The  location  is  at  depth  z  =  1  [mm]  below  where  the  static  solution  has  a  maximum  (see  Fig.  3).  As 
in  the  case  of  the  mean  normal  pressure,  we  observe  that  the  solutions  are  identical  up  to  a  critical 
frequency,  given  by  ucr  =  kcr/L  or  while  /  < 


standard  linear  solid.0  The  transfer  function  for  /(  is  for  this  model, 


M(s)  =  no 


1  4“  trs 
1  +  ^ trS 


where  no  is  the  steady  state  response,  tr  the  retardation  time  and  a  is  the  overshoot  of  a  step 
response.  A  schematic  representation  of  the  model  as  dampers  and  springs  is  shown  in  Fig. 
4-5a,  and  the  step  response  in  Fig.  4-5b,  i.e.  the  response  in  stress  to  a  step  in  strain. 


This  model  can  be  generalized  to  model  materials  with  more  than  one  time  constant.  The 
transfer  function  for  the  generalized  model  is 


n 

m0) = mo  n 

X— 1 


1  ~f~  S 

1  +  ^tTis' 


sAlso  called  Kelvin  model. 


(a)  3-Parameter  Standard  Linear  Solid 


(b)  Response  of  stress  to  a  step  in  strain 


Figure  4-5:  (a)  Shows  a  diagram  of  the  3-parameter  linear  standard  solid  model  using  spring  and  damper 
elements,  (b)  Shows  the  response  in  stress  to  a  step  change  in  strain  for  the  3-parameter  linear  solid 
defined  in  Fig.  4-5.  The  relation  between  the  parameters  is  that  p0  =  &i k2/(ki  +  k2),  a  =  (A,’i  +  k-2)/k-2, 
and  tr  —  b\/(ki  k2). 


The  transfer  function  approach  to  viscoelasticity  provides  a  powerful  tool  for  the  modeling 
of  the  viscoelastic  properties  of  different  materials.  One  still  needs,  however,  to  estimate  the 
parameters  in  the  transfer  function  for  each  material  and  under  conditions  that  are  as  close  as 
possible  to  the  conditions  under  which  the  model  is  to  be  applied. 


This  model  captures  the  general  behavior  of  an  viscoelastic  solid,  in  that  the  net  effect  is  a 
phase  delay  between  the  loading  and  the  displacement  occuring  because  energy  is  absorbed  by 
the  medium.  This  phase  delay  is  one  of  the  factors  that  makes  a  dynamic  formulation  necessary, 
even  when  only  viscoelastic  effects  are  dominant  and  the  encoding  solution  is  identical,  in  form, 
to  the  static  solution.  This  we  will  discuss  in  the  section  3. 


3.  The  Decoding  Problem 


Tactile  information  is  obtained  with  mechano-sensors  embedded  in  a  deformable  medium  that 
is  in  contact  with  a  shaped  object.  Mechanosensors,  embedded  within  the  medium  can,  at 
best,  provide  information  on  skin  surface  shape  and  surface  traction  distribution.  From  this 
information,  contact  region,  shape  of  the  object  and  the  contact  state  need  to  be  inferred. 
The  problem  at  hand  is  therefore  the  decoding,  i.e.  reconstruction  of  the  surface  shape  of 
the  medium,  the  tractions  on  the  surface  as  well  as  the  region  of  contact  using  sub-surface 
mechanosensory  information. 
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An  important  part  of  the  decoding  solution  is  the  selection  of  the  sensor  signals,  i.e.  the 
particular  strains  or  stresses  the  sensors  should  pick  up.  As  the  introduction  of  the  dynamics 
does  not  alter  the  symmetry  properties  of  the  halfspace  we  will  refer  to  the  treatment  of  the 
static  solution  in  chapter  3  and  not  need  repeat  it  here. 

In  Eq.  (4.4),  we  have  chosen  to  represent  the  solution  for  the  encoding  problem  in  terms  of 
transfer  functions  in  both  the  spatial  and  temporal  frequency  domains.  We  will  show  in  this 
section  that  for  the  solution  to  the  decoding  problem,  it  is  more  efficient  to  use  state-space 
representation  for  the  temporal  dimension  when  possible. 


3.1  The  Decoding  Solution 


The  solution  to  the  decoding  problem  essentially  gives  us  an  estimate  of  the  conditions  on  the 
surface  given  a  time  history  of  the  sensor  signals.  This  estimate  must  be  such  that  it  minimizes 
some  error  measure.  One  such  measure  is  the  expected  total  square  error  defined  as 


£  =  E 

=  E 


f  r 

/  /  IlCOc.JbO  -  Z(x,y,t)W2dxdy 

J  J  —  oo 

J J  t)  ((Ai  Ct/'y  ,  t  )  ||  du!XdU/y 


(4.7) 


where  £  and  £  are,  respectively,  the  true  and  estimated  signals  describing  the  conditions  at  the 
surface.  The  latter  relation  follows  from  Parceval’s  theorem,  and  it  shows  us  that  we  can  with 
identical  results  minimize  the  error  in  the  domain  of  spatial  frequency. 


When  the  effects  of  inertia  is  included  in  the  encoding  solution,  the  resulting  TFMs  are 
irrational  in  the  frequency  variable  s.  This  indicates  that  the  system  being  described  does  not 
have  a  finite  number  of  states  and  hence  cannot  be  represented  using  the  state-space  approach. 
We  have  shown  in  section  2.2  that  the  effects  of  inertia  are  only  of  importance  for  higher 
temporal  frequencies  which  are  in  most  cases  out  of  the  range  of  interest  in  most  applications. 
Hence,  we  will  concentrate  on  the  solution  when  viscoelastic  effects  dominate. 


The  decoding  problem  as  posed  here  is  in  most  aspects  identical  to  optimal  linear  filtering 
problems  encountered  in  signal  processing  and  control  theory.  In  those  fields,  extensive  research 
has  lead  to  the  now  accepted  conclusion  that  for  this  problem  effective  solutions  can  only  be 
obtained  for  signals  that  can  be  described  using  rational  power  spectral  functions  [86,  37,  52]. 
The  system  can  then  be  represented  in  a  state-space  and  the  resulting  filter  is  the  wellknown 
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Kalman  filter. 


Translating  the  shear  modulus  fi(s)  from  transfer  function  to  state  space  representation  can 
be  done  using  any  of  the  standard  control  theory  methods,  cf.  [62].  In  order  to  represent  the 
surface  load  to  sensor  TFM ,(TS/  in  Eq.  (3.4)),  in  state  space,  we  need  to  use  as  many  state 
variable  as  the  order  of  p(s)  dictates  for  each  of  exz,  eyz.  However,  since  pn  is  a  stress  variable 
and  the  input  is  load,  p,  does  not  appear  in  the  corresponding  TFs.  Therefore  for  pn,  there  are 
no  viscoelastice  effects  that  need  to  expressed  as  state  variables.  Effectively,  we  assume  that 
the  response,  of  the  mean  stress  sensors  to  surface  load,  is  instantaneous.  Denoting  the  states 
for  exz  and  eyz  as  ££xz,  ££yz  €  IR”,  respectively,  where  ??,  is  the  order  of  p(s),  we  get,  using 
control  canonical  form  for  each  sensor  variable,  that  the  TFM  can  be  represented  in  state  space 
as 
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where  Ay  6  IR"xn  and  c[L  €  IR"  are  obtained  from  p(s),  and  the  input  matrix  is  made  of  the 
static  encoding  solution  TFs  from  Tsf  in  Eq.  (3.4)  ’as  follows 
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The  Kalman  Filter 


In  order  to  formulate  the  decoding  problem  as  a  filtering  problem  we  assume  that  known 
properties  of  the  surface  signal  can  be  modeled  using  a  colored  noise  signal,  i.e.  a  signal  filtered 
through  a  known  linear  system  that  characterizes  the  expected  spatial  and  temporal  properties 
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of  the  signal,  primarily  bandwidth.  The  form  and  the  parameters  of  that  system  need  to  be 
assessed  for  each  model  or  implementation,  but  this  paradigm  provides  a  tool  to  incorporate 
and  utilize  known  properties  of  the  input  in  the  decoding.  For  example,  if  we  assume  that  the 
temporal  and  spatial  characteristics  are  separable  and  write 

E[f{ujx,LJy,s)f  (cnx,u;y,s)]  =  Ps(wx, wy)$(s), 

where  Ps  is  the  spatial  power  spectrum  matrix  of  the  solution  discussed  in  chapter  3  section 
3.3  but  $  is  the  temporal  power  spectrum  matrix  that  can  be  used  as  a  shaping  filter  to  model 
the  temporal  characteristics  of  the  load  signals.  $  can  be  any  valid  rational  power  spectrum, 
but  a  simple  example  of  such  a  model  is  the  Markov  process  model,  where  the  power  spectrum 
and  the  corresponding  shaping  filter  is  given  by 

^ ,  w  JTT  vw 

$00  =  — TTln  WuadH*  =  ——5 
—  sz  +  p*  s  +  p 

where  0  represents  the  temporal  bandwidth  of  the  load  signal.  Using  the  above  example,  the 
spatial  and  temporal  spectrum  characteristics  can  then  be  combined,  e.g.  for  the  normal  load 

7, 

7z  =  -P7  z  +  \j20Pszz(ujx,ujy)tu(t) 

where  w  is  a  white  noise  signals,  Pzz  is  the  spatial  power  spectrum  of  fz,  and  we  have  made 
the  simplifying  assumption  that  Ps  is  diagonal. 


Having  expressed  both  the  encoding  solution  and  the  known  spectral  properties  as  state 
variables,  we  combine  those  and  obtain  the  system  in  form  that  we  can  directly  apply  the 
Kalman  filter  to 


£  =  A£  +  Bw 
C  -  Cl  +  v. 


(4.9) 


.4,  B  and  C  are  obtained  by  augmenting  the  system  Eq.  (4.8)  with  the  states  used  to  model 
the  temporal  spectral  properties  of  the  surface  loads  or  or  displacements  as  discussed  above. 
The  resulting  state  is  £  and  we  have  used  £  for  the  measurement  vector,  i.e.  £  =  [pn,  exz,  eyz]T 
The  input  iu  and  the  sensor  noise  v  are  assumed  to  be  white  noise  processes  with  the  following 
properties: 

E[w(t)w(r)*]  =  Q6(t-r) 

E[v{t)v{r)*}  =  R6(t  —  t)  (4.10) 

E  [uJ(t)?(r)*]  =  0. 

where  R  =  Pn  the  same  as  delineated  in  section  3.3,  chapter  3  and  Q  =  Ps,  if  temporal  and 
spatial  properties  are  separable,  as  discussed  above.  We  further  note  that  the  Kalman  filter 
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allows  A,B,C,Q  and  R  functions  of  spatial  frequency  and  time,  i.e.  A  =  A(u>x,uy,t),  hence 
the  separability  assumption  above  is  not  necessary  but  we  only  present  it  to  show  the  relation 
of  this  solution  to  the  static  solution  (see  also  Appendix  E). 


For  linear  systems  in  state  space  it  is  known  that  a  filter  that  minimizes  the  expected 
covariance  of  the  estimation  error  P  —  E  tr  ((£  —  £)(£  —  0*)]  and  hence,  minimizes  the  error 
in  Eq.  (4.7),  is  the  Kalman  filter  given  by  the  following  equations: 


?  =  a}  +  k(c-co 

K  =  PC*R~l  (4-n) 

P  =  AP  +  PA*  -PC'R-'CP  +  BQB*  P(0)  =  P0 

and  we  get  our  solution  to  the  decoding  problem  from  the  state  estimate  since  the  inputs  are 
contained  in  that  part  state  vector  that  models  the  input  characteristics.  (For  examples  see 
simulations  in  section  4.. 


This  state  representation  in  the  spatial  frequency  domain  is  a  straightforward  and  a  powerful 
extension  to  the  static  decoding  solution.  It  retains  the  numerical  efficiency  of  the  Fourier 
transform  (i.e.  with  the  use  of  FFT)  and  at  the  same  time  provides  a  well  developed  modeling 
tool.  The  Kalman  filter  further  has  many  guaranteed  properties  and  an  important  one  is  that 
if  the  signals  involved  (w  and  v)  can  be  assumed  to  be  Gaussian  the  the  Kalman  filter  is  the 
best  possible  filter  for  a  large  class  of  criterion  [6],  i.e.  it  can  be  shown  that  it  minimizes  the 
expectation  of  nondecreasing  functions  of  the  magnitude  of  the  error. 


Remark:  Many  of  the  assumptions  made  above  (e.g.  that  E  | iu(t)vT(r)  =  0)  can  be  dropped 
or  weakened  but  we  believe  that  more  detailed  discussion  on  the  properties  of  the  Kalman  filter 
is  beyond  the  scope  of  this  thesis  and  the  reader  is  refered  to  the  rich  literature  that  exist  on 
the  kalman  filtering  and  its  many  applications,  (e.g.  [6,  37].) 


Remark:  Even  though  this  formulation  of  a  solution  to  the  decoding  problem  appears  quite 
different  than  the  one  obtained  in  the  case  of  static  loading  in  chapter  3  it  is  in  fact  equivalent 
in  the  sense  that  the  Kalman  filter  solution  can  be  shown  to  converge  to  the  Wiener  filter  under 
static  conditions.  This  we  show  in  appendix  E. 
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4.  Simulations 


The  simulations  presented  here  are  designed  to  clearly  draw  out  the  specific  aspects  of  the  en¬ 
coding  and  decoding  problems  treated  in  this  chapter,  i.e.  the  dynamic  effects  of  viscoelasticity 
on  these  problems.  In  order  not  to  clutter  the  results,  we  chose  the  examples  to  be  as  simple  as 
possible  while  still  being  able  to  reach  the  main  objective  of  this  simulation  section,  which  is  to 
confer  the  following:  1)  Viscoelastic  effects  need  be  accounted  for  if  one  wants  to  correctly  solve 
the  encoding  and  decoding  problems  under  dynamic  conditions.  2)  State  space  representation 
and  the  Kalman  filter  can  be  successfully  applied  to  solve  the  above  problems. 


The  example  presented  here  will  be  a  model  of  a  smooth  rigid  body,  shaped  in  the  form  of  a 
paraboloid,  indenting  a  viscoelastic  halfspace.  The  normal  load  is  assumed  to  be  a  time- varying 
Hertzian  distribution,  which  is  one  of  the  few  available  analytical  solutions  to  viscoelastic  con¬ 
tact  problems.  The  time  varying  Hertz  distribution  can  be  shown  to  correctly  model  the  in¬ 
dentation  of  an  object  of  this  shape  into  a  viscoelastic  halfspace  during  the  indentation  phase, 
i.e.  while  the  contact  area  is  nondecreasing  [44,  30]. 


The  contact  area  is  time-varying  but  the  contact  radius  is  obtained  by  integrating  the 
relation 


3  RFz(s) 
8  fi(s) 


(4.12) 


where  the  parameter  R  is  the  curvature  of  the  object,  as  u*  =  (x2  +  y2)j'2R  and  we  have 
assumed  that  the  material  is  incompressible  {v  =  0.5).  We  use  here  the  simple  Kelvin  model, 
discussed  in  section  2.3  for  the  shear  modulus,  or 


h(s) 


1  +  5 

1  +  s/3 


i.e.  tr  =  1  s  and  ar  =  3.  We  further  choose  the  total  normal  load  trajectory  to  be  a  low-pass 
filtered  step  (see  Fig.  4-6a)  The  resulting  time  history  for  the  contact  radius  a(t)  is  shown  in 
Fig.  4-6b  and  corresponding  surface  load  and  displacement  distributions  are  shown  in  Figs. 
4-6c  and  4-6d  for  the  time  points  t  =  0.02,  0.1,  1,  and  4  s.  The  surface  normal  load  and  the 
displacement  distributions  shown  in  Fig.  4-6  are  calculated  using  the  known  Hertzian  solution 
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Figure  4-6:  (a)  Shows  the  time  history  of  the  total  normal  load  Fz(t)  which  is  a  low-pass  filtered  step, 
(b)  Shows  the  corresponding  contact  radius  time  history  obtained  from  Eq.  (4.12).  (c)  Shows  sections 
of  the  resulting  surface  normal  load  distribution,  /2, obtained  using  Eq.  (4.13)  at  for  different  times, 
(t  —  [0.02,  0.1,  1,  4]).  We  note  that  the  maximum  load,  is  not  monotonic  in  time  as  the  total  load  but 
reaches  a  maximum  in  t  €  [0.1, 1]  (see  also  Fig.  4-7.  (d)  Shows  sections  of  the  resulting  surface  normal 
displacement  for  the  same  time  points  as  (c). 
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Although  the  solution  presented  in  section  3.  is  valid  for  a  general  3D  loading  we  only  include  the 
normal  load  since  we  are  concentrating  on  the  dynamic  and  viscoelastic  effects  in  this  chapter 
and  they  are  independent  of  the  effects  of  the  shear  loads.  Hence,  the  results  and  comments  on 
the  effects  of  static  shear  loading  obtained  in  chapter  3  are  equally  applicable  when  dynamics 
and  viscoelasticity  are  included  and  we  therefore  do  not  repeat  them  here. 


4.1  The  Encoding 


We  pose  the  encoding  problem  here,  as  the  prediction  or  calculation  of  the  subsurface  signals 
that  the  sensors  transduce,  given  the  surface  load  or  displacement  distributions.  We  have 
already  decided  that  the  sensors  are  the  set  (pn,  exz,  eyz)  and  since  we  later  on  intend  to  apply 
the  Kalman  filter,  we  use  the  state  space  representation  given  in  Eq.  (4.8).  The  state  space 
matrices,  not  already  specified,  for  this  system  become 

=  and  Cu  =  tr  (l  — —  ) 

tr  \  Clr  / 

and  Fig.  4-7  shows  the  shape  of  the  resulting  signals  for  the  same  time  points  as  Fig.  4-6.  We 
note  that  the  shape  of  mean  normal  stress,  pn  sections  (Fig.  4-7a-b)  is  a  low-pass  filter  version 
of  the  surface  normal  load  in  Fig.  4-6c  and  that  the  response  is  instantaneous.  The  shear  strain 
response  shown  in  Fig.  4-7c,  on  the  other  hand,  is  not  instantaneous  and  there  the  viscoelastic 
effects  become  clearly  visible.  The  instantaneous,  or  static,  response  is  shown  as  a  dotted  line 
and  we  see  that  it  behaves  quite  differently,  with  time,  than  the  viscoelastic  response,  shown  as 
a  solid  line.  We  further  note  that  the  sections  corresponding  to  t  =  4  s  coincide  so  the  responses 
converge  in  the  end  as  expected. 
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Figure  4-7:  (a)  Shows  the  time  history  of  maximum  mean  normal  pressure,  max(pn).  The  time  points 
selected  to  show  the  distribution  profiles  are  marked  with  x  on  the  plot  and  we  note  that  there  is  a 
maximum  at  t  ss  0.3s.  (b)  Shows  sections  of  the  mean  normal  pressure  distribution,  (c)  Shows  sections 
of  the  zs-shear  strain  distribution.  The  solid  lines  show  the  viscoelastic  response  but  the  dotted  lines 
what  the  predicted  response  would  be  if  viscoelastic  effects  were  ignored.  We  note  that  they  are  quite 
different  for  all  time  points,  except  for  t  =  4  s  when  the  response  profiles  coincide  as  then  the  static 
assumption  becomes  valid. 
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4.2  The  Decoding 


Here  the  problem  is  to  reconstruct  the  surface  load  and  displacement  distributions  from  the 
sensor  signals.  We  only  consider  normal  loads  so  the  input  signal  is  ID  and  we  therefore  only 
need  one  sensor  signal  for  which  we  choose  pn.  We  could  use  the  other  signals,  i.e.  exz  and 
eyz  instead  or  as  well  but  that  would  not  have  any  effects  on  the  dynamic  side  of  this  problem, 
which  we  want  to  expose. 

Reconstructing  the  surface  load  distribution,  /*,  fromp„  will  not  involve  viscoelastic  effects 
as  our  assumptions  implicitly  say  that  this  response  is  instantaneous.  The  reconstruction  of  the 
displacement  distribution,  uz,  does  however  involve  dynamics  and  then  we  use  the  above  state 
space  model  with  the  transfer  function  TpnU.  =  2/ufiexp(— Qz)  in  both  the  B  and  C  matrices. 
We  further  use  the  Markov  process  to  model  the  temporal  properties  of  the  surface  signals 
and  choose  the  parameter  0  to  be  large  enough  to  match  approximately  the  bandwidth  of  the 
surface  signals.  Finally,  we  add  a  spatially  high-pass  filtered  Gaussian  sensor  noise,  so  that  the 
rms.  error  to  signal  ratio6  is  ss  23%,  see  Fig.  4-8. 

The  results  from  the  Kalman  filter  are  shown  in  Figs.  4-9a-c  and  there  we  note  that  the 
filter  successfully  estimates  the  surface  signals  and  that  the  rms.  error  ratios  are  ~  14%  for  /-, 
the  surface  normal  load  and  3%  for  uz,  the  surface  normal  displacent. 

In  this  section  we  have  given  an  example  of  how  the  theory  presented  in  this  chapter  can  be 
applied  to  solve  the  encoding  and  decoding  problems  under  dynamic  conditions.  The  results  of 
the  simulations  indicate  that  the  state  space  approach  can  be  used  to  model  viscoelastic  prop¬ 
erties  and,  simultaneously  include  known  spatial  and  temporal  properties  of  the  signals.  The 
encoding  results  indicate  that  viscoelastic  effects  can  have  significant  effects  on  the  predicted 
sensor  response  and  that  viscoelastic  effects  need  to  be  accounted  for  under  dynamic  conditions. 
The  decoding  results  give  us  that  the  Kalman  filter  is  successful  in  filtering  out  sensor  errors 
and  gives  relatively  accurate  estimates  of  the  surface  signals,  where  the  displacement  estimates 
are  more  accurate  than  the  load  distribution  estimates. 

6That  is,  square  root  of  the  ratio  between  the  integrated  squared  error  and  squared  signal,  respectively 
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(a) 


(b) 


Figure  4-8:  (a)  Shows  the  distribution  of  the  mean  normal  stress,  pn,  at  sensor  depth,  (b)  Shows  the 
same  distribution  with  spatially  high-pass  filtered  measurment  noise  added.  The  rms.  error  level  is 
«23%. 
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Figure  4-9:  The  results  of  decoding  simulations,  (a)  Shows  the  actual  and  decoded  estimates  of 
surface  normal  displacement  distribution  uz.  The  estimates  are  shown  using  dotted  lines  but  the  actual 
displacements  using  solid  lines,  (b)  Shows  the  actual  and  decoded  estimates  of  surface  normal  load 
distribution  fz.  The  estimates  are  shown  using  dotted  lines  but  the  actual  loads  using  solid  lines,  (c) 
Shows  the  rms.  errors  as  percentage  ratios  to  the  total  response  of  each  signal.  We  note  that  both 
the  load  and  displacement  errors  are  less  than  the  sensor  noise  and  that  the  displacment  error  ratio  are 
always  less  than  the  load  error  ratio  but  that  is  caused  by  the  smoothing  nature  of  continuous  materials. 
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5.  Conclusions 


In  this  chapter  we  have  delineated  the  aspects  of  haptic  computational  theory  that  relate  to 
dynamic  properties  of  the  fingerpad.  Separating  the  problem  of  shape  identification  into  an 
encoding  and  decoding  problem  we  give  a  solution  to  the  encoding  problem  using  a  Transfer 
function  Matrix  formulation  obtained  using  a  combined  2D-Fourier-Laplace  transform.  We  show 
how  the  effects  of  inertia  and  viscoelasticity  enter  the  solution  and  analyze  the  contribution  of 
each.  We  further  show  that  the  effects  of  inertia  can  be  neglected  for  a  significant  lower  range 
of  the  temporal  (lower  bandwidth  signals).  We  give  a  solution  to  the  decoding  problem  in 
the  form  of  a  Kalman  filter  operating  on  the  2D-Fourier  transform  of  the  sensor  signals.  The 
Kalman  filter  is  a  numerically  effective  and  well  developed  tool  with  many  desirable  properties, 
and  it  has  proven  its  validity  in  numerous  applications.  Finally  we  show  how  the  solutions  to 
the  encoding  and  the  decoding  problems  become  identical  to  the  solutions  to  the  static  problem 
given  in  chapter  3  when  the  surface  loads  become  static  or  quasi-static. 

A  computational  model  of  dynamic  tactile  contact  provides  an  important  extension  to  pre¬ 
vious  results  on  the  static  haptic  identification  problem.  Robots  can  generally  be  expected  to 
operate  under  dynamic  conditions  where  contact  loads  cannot  be  assumed  to  be  quasi-static.  It 
is  also  known  that  for  humans  stroking  enhances  performance  when  carrying  out  haptic  tasks. 
We  further  know  that  the  characteristic  time  constants  of  the  fingerpad  are  on  the  order  of  a 
second  or  more  and  that  viscoelastic  effects  are  significant  and  cannot  be  ignored.  Finally  in 
virtual  environments,  it  is  the  dynamic  interaction  that  simultaneously  provides  most  of  the 
illusion  and  implementation  difficulties,  making  detailed  modeling  necessary.  A  related  problem 
is  how  the  information  that  is  obtained  using  the  methods  described  in  this  chapter  can  be  used 
to  touch  and  grasp,  and  consequently  explore  environments  or  manipulate  objects  ,  but  that 
will  be  the  subject  of  next  chapter. 
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Chapter  5 


Exploration  Strategies  and  Finger 
Control 


1.  Introduction 


A  comprehension  of  the  nature  of  contact,  grasping  and  subsequent  exploration  or  manipulation 
is  of  significant  importance  in  both  physiological  and  robotics  research.  In  physiological  research 
the  objective  is  to  understand  the  mechanics  of  hand  movements,  especially  which  factors  govern 
the  shape  of  motion  trajectories  and  more  importantly  why  they  do  so.  This  knowledge  has  in 
addition  to  its  scientific  value,  uses  in  such  areas  as  rehabilitation,  learning  and  development  of 
motion  skills  and  ergonomic  studies,  just  to  name  a  few.  On  the  robotics  side  the  goal  has  been 
to  make  a  robots  perform  tasks  that  demand  human  level  of  dexterity.  This  goal  has  proven 
elusive,  mainly  because  of  the  inherent  complexity  of  the  mechanics  involved.  But  the  difficulty 
is  easy  to  underestimate  considering  the  relative  ease  humans  seem  to  have  with  performing 
these  tasks.  Robots  have  until  now  also  lagged  behind  their  human  counterparts  on  the  sensory 
side  and  it  is  only  recently  that  sensor  arrays  have  been  produced  that  in  some  ways  duplicate 
the  properties  of  the  human  sensory  organs.  The  problem  of  grasping  and  manipulating  an 
object  has  been  investigated  by  many  researchers  (for  review  see  [69])  and  in  view  of  that  and 
the  above  we  will  in  this  chapter  discuss  how  this  sensory  information  can  be  used  for  the 
exploration  of  objects.  We  will  try  to  make  our  arguments  independent  of  implementation  and 
therefore  equally  applicable  to  robots  and  humans.  By  this,  we  hope  to  extract  a  common 
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thread  in  a  web  that  could  be  called  a  computational  theory  of  haptics. 


In  chapters  3  and  4  we  have  discussed  the  processing  of  sensor  signals  resulting  in  algo¬ 
rithms  and  guidelines  for  sensor  design  that  together  give  optimal  estimates  of  the  loading  and 
displacement  distributions  on  the  surface  of  the  fingerpad.  That  is  the  Surface  Signal  Identifi¬ 
cation  Problem ,  which  details  can  be  found  in  chapters  3  and  4.  Fig.  5-1  shows  a  diagram  which 
explains  the  relation  of  that  part  to  the  problems  we  are  going  to  discuss  here.  In  this  chapter 
we  will  focus  on  two  blocks  of  Fig.  5-1,  the  Exploration  Strategy  and  the  Finger  Control  blocks. 
Our  main  goal  is  to  explore  objects,  but  the  Exploration  Strategy  block  refers  to  the  problem 
of  inferring  object  properties  such  as  shape,  texture,  and  compliance  and  the  state  of  contact 
from  the  estimated  surface  signals.  This  involves  determining,  given  a  property  that  is  to  be 
explored,  what  kind  of  sensor  information  and  also  what  kind  of  action  is  needed  for  successful 
exploration.  Thus,  in  addition  to  actual  information  processing,  a  specific  action,  that  brings 
the  desired  information  to  the  sensors,  is  also  needed.  This  is  the  role  of  the  Finger  Control 
block,  but  it  refers  to  the  translation  of  the  desired  actions  into  command  trajectories  for  the 
fingerpad.  We  will  show’  that  for  continued  exploration  we  will  need  information  about  the 
local  curvature  of  the  object  at  the  contact  location.  We  further  give  a  method  to  extract  that 
information  from  the  surface  signals  estimates  in  a  robust  manner.  Finally,  Fig.  5-1  shows  how 
the  command  trajectories  are  fed  into  a  low-level  manipulator  controller  which  we  assume  is  ca¬ 
pable  of  producing  motor  commands  that  will  make  the  manipulator  move  accurately  enough1 
along  the  command  trajectories,  which  in  turn  closes  the  loop  at  the  fingerpad. 

The  Supervisor  figure  in  Fig.  5-1  provides  extermal  information  labeled  control  strategy 
which  for  example  includes  the  desired  direction  of  stroking.  We  therefore,  assign  to  the  Super¬ 
visor  the  role  of  providing  ‘global’  information  which  is  not  readily  accessable  through  tactile 
channels,  e.g.  visual  information.  One  of  the  implications  of  this  assumption  is  that  we  will 
not  dwell  on  various  problems  that  arise  when  an  object  is  explored,  such  as  how  can  we  ensure 
we  will  cover  the  whole  surface  of  an  object  when  stroking  it.  We  believe  that  such  questions 
are  more  naturally  answered  using  either  heuristic  or  AI  techniques  and  as  such  are  beyond  the 
scope  of  this  thesis. 

We  are  not  aware  of  any  reference  that  treats  the  problems  discussed  here  as  whole  and 
we  believe  the  reason  is  the  fact  that  it  has  been  addressed  from  different  directions  in  various 
disciplines.  Exploration  strategies  have  mostly  been  investigated  from  the  human  side.  [42,  43] 


1The  manipulator  we  discuss  here  can  be  a  human  arm  or  finger  as  well  as  a  robot  manipulator.  We  acknowl¬ 
edge  that  the  human  arm  has  in  general  inferior  position  and  force  accuracy  when  compared  to  robot  arms. 
However;  we  also  know  that  measured  in  dexterity  it  is  superior  so  the  human  accuracy,  or  lack  thereof,  must  be 
sufficient  in  some  way  and  this  is  what  we  mean  by  saying  “accurately  enough”. 
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investigates  exploratory  procedures,  [48]  delineates  sensor  channels  and  coding  in  general,  but 
[76,  77,  48, 12,  33,  29]  discuss  the  same  with  respect  to  specific  object  properties,  mainly  texture. 
Most  of  the  definitions  in  these  references  are  implementation  independent  and  therefore  usable 
for  robots  as  well.  It  is  only  the  definition  of  object  properties  such  as  texture  that  is  often 
subjective  and  often  hard  to  get  any  consensus  on  what  the  term  texture  means.  We  will 
therefore  explicitly  define  the  object  properties  we  will  discuss  so  that  they  can  be  quantitively 
identified.  For  shape  sensing  we  rely  mainly  on  image  analysis  literature,  [1,  31,  3,  50,  78,  27], 
and  for  texture  measures  w'e  borrow  concepts  from  tribolobj'  in  work  on  surface  roughness 
[79,  47,  46,  58,  25].  The  robotics  literature  dominates  in  the  formulation  of  control  laws  as  can 
be  expected.  For  contact  and  grasping  stability  during  manipulation  we  confer  [69,  54,  7],  but 
for  stroking  and  manipulation  kinematics  we  confer  [57,  53,  55],  and  finally  for  the  formulation 
of  control  laws  we  confer  [2,  26,  11]. 

The  contributions  of  this  chapter  lie  mainly  in  the  arrangement  in  a  common  framework 
the  different  components  and  processes  of  haptic  exploration.  These  components  and  processes 
have  until  now  been  treated  treated  separately,  but  we  identify  them  as  different  facets  of  a 
larger  structure,  ie.  the  structure  illustrated  in  Fig.  5-1.  On  the  more  detailed  level,  we  define 
and  discuss  in  a  consistent  manner  the  constituting  parts,  i.e.  the  Exploratory  Procedures  for 
each  object  Property  and  contact  state  and  also  the  relevant  Sensor  Channels  and  Codes  in 
each  case.  We  finally  show  how  these  exploration  strategies  can  be  implemented  and  specifically 
delineate  how  stroking  can  be  achieved  using  robust  estimates  of  local  shape  information. 

This  chapter  is  organized  sectionwise  as  follows:  In  section  2  we  discuss  the  Exploration 
Strategies  including  definitions  of  the  term  used,  i.e.  Exploratory  Procedures,  Sensor  Channels 
and  Sensor  Codes.  We  then  apply  the  tools  introduced  for  the  identification  of  object  properties, 
and  the  state  of  contact.  In  section  3  we  discuss  translation  of  the  Exploratory  Procedures  into 
control  laws  and  we  do  that  separately  for  directions  tangential  and  normal  to  the  fingerpad. 
We  finally  give  our  conclusions  and  remarks  in  section  4. 


2.  Exploration  Strategies 


In  this  section  we  will  discuss  exploration  strategies  suitable  for  object  property  or  state-of- 
contact  identification.  More  specifically  we  will  list  the  types  of  actions,  i.e.  Exploratory 
Procedures  human  and  robots  have  available.  We  will  also  list  different  types  of  sensors,  ie. 
Sensor  Channels  and  further  what  part  of  the  sensor  signal  is  of  interest,  i.e.  Sensor  Codings. 
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Table  5.1:  Explicitly  lists  the  available  exploratory  procedures,  sensor  channels  and  codes  as  well  as 
object  properties.  For  further  clarifications  of  terms  used  in  the  table  see  text. 


Exploratory 

Procedures 

Sensor 

Channels 

Sensor 

Codes 

Static  Touch 

SA1 

Intensive 

Indenting 

Sensor  Array 

Temporal 

Stretch 

RA1 

Spatial 

Stroking 

SA2 

Spatio- 

Tapping 

PC 

temporal 

(Enclosure) 

Accelerometer 

Object 

Shape,  Texture,  //Texture 

Properties 

Compliance,  State  of  Contact 

Finally,  we  will  list  and  define  the  object  properties  that  are  identifiable  using  the  previously 
listed  exploratory  procedures  and  sensor  channels  and  codes.  For  overview,  Table  5.1  explicitly 
lists  the  available  exploratory  procedures,  sensor  channels  and  codes  as  well  as  the  object 
properties  that  we  discuss  in  more  detail  below. 


Exploratory  Procedures 


The  term  Exploratory  Procedure  as  defined  in  [42],  refers  to  a  “stereotyped  movement  pattern 
having  certain  characteristics  that  are  invariant  and  others  that  are  highly  typical".  Here  we 
will  only  include  simple  movements  that  only  require  the  movement  of  a  single  fingerpad  and 
exclude  more  complicated  procedures  such  as  function  tests  or  part  motion  tests  that  may 
require  the  relative  motion  of  two  hands.  The  “enclosure”  procedure  actually  violates  this 
condition  since  it  would  require  multiple  fingerpads,  but  we  include  it  since  it  adds  insight 
when  explaining  object  shape  identification  (see  section  2.1). 


1.  Static  Touch  is  the  most  basic  procedure  and  it  can  only  give  information  on  the  shape 
of  the  fingerpad  within  the  contact  area,  and  only  spatial  information  is  available. 

2.  Indenting  gives  also  information  about  the  force-displacement  relationship  of  an  object 
and  can  therefore  be  used  to  infer  the  normal  compliance  of  the  object. 
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3.  Stretch  gives  same  addition  as  indenting  but  in  the  tangential  direction  and  will  therefore 
give  information  on  the  lateral  stiffness  of  an  object. 

4.  Stroking  adds  the  temporal  dimension  and  the  spatio-temporal  information  becomes 
available.  One  feature  of  stroking  is  that  if  position  accuracy  is  better  than  sensor  resolu¬ 
tion  then  one  can  increase  the  effective  spatial  bandwidth  as  then  the  surface  is  effectively 
re-sampled.  The  extents  of  this  bandwidth  enhancement  are  limited  by  the  position  ac¬ 
curacy. 

5.  Tapping  is  a  ballistic  movement  where  the  finger  impacts  the  object,  and  provides  the 
same  information  as  indenting. 

6.  Enclosure  is  the  simultaneous  contact  of  multiple  fingers  and  may  include  other  parts 
of  the  hands  such  as  the  palms. 


Physiological  experiments  involving  analysis  of  videotaped  hand  movements  indicate  that  hu¬ 
mans  use  these  exploratory  procedures  to  a  high  degree  and  further  that  they  select  different 
procedure  to  match  the  required  object  property  [43].  We  therefore  believe  that  this  approach 
is  of  value  and  can  be  used  for  robotic  exploration. 


Sensor  Channels  and  Codes 


The  concept  of  sensor  channels  and  the  coding  principle  introduced  in  [34,  40]  has  its  roots  in 
neurophysiology  where  researchers  were  interested  in  how  different  textures  and  other  stimuli 
to  the  fingerpad  of  monkeys  was  coded  in  the  measured  peripheral  neural  signals.  There  this 
classification  has  been  proven  to  be  useful  in  generating  hypothesis  on  phenomena  such  as 
roughness  perception.  It  can  however,  also  be  useful  for  the  identification  of  other  object 
properties  such  as  shape  and  compliance. 


Sensor  Channels 


The  need  to  define  different  sensor  channels  arises  from  the  fact  that  there  are  known  at  least 
four  different  types  of  mechanoreceptors  the  human  fingerpad.  These  are  frequently  classified 
as:  slowly  adapting  type  I  (SA1),  slowly  adapting  type  II  (SA2),  rapidly  adapting  (RAl)  and 
Pacinian  (PC).  We  have  in  table  5.1  where  applicable  identified  the  robot  sensor  types  that  are 
roughly  parallel  to  the  human  mechanoreceptors.  The  criteria  we  used  in  matching  human  and 
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Table  5.2:  Lists  the  different  sensor  channels  for  humans  and  robots  and  matches  them  according  to 
temporal  and  spatial  bandwidth. 


Human 

Meclianorepectors 

Spatial  and  Temporal 
bandwidths 

Robotic 

Sensors 

SA1 

Low  temporal,  high  spatial 

Strain  gauge  array 

SA2 

Low  temporal,  low  spatial 

Force  sensor 

RAl 

Medium  temporal,  medium  spatial 

PC 

High  temporal,  low  spatial 

Force  sensor 

robotic  are  the  sensor’s  differing  temporal  and  spatial  bandwidth  of  the  sensors.  Table  5.2  lists 
how  this  comparison  is  made. 

Sensors  with  limited  spatial  and  temporal  frequency  range  are  easier  to  make  than  sensors 
that  cover  large  ranges  of  both.  It  is,  however,  possible  to  make  several  sensors  each  targeting 
a  different  frequency  range  and  combine  the  information  from  the  different  sensor  channels  to 
cover  the  frequency  ranges  of  interest.  How  this  concept  is  realized  in  the  human  fingerpad 
is  shown  schematically  in  Fig.  5-2,  where  the  the  temporal  and  spatial  frequencies  are  on  the 
horizontal  and  vertical  axis,  respectively. 


An  interesting  coupling  between  the  temporal  and  spatial  frequency  responses  occurs  during 
stroking  and  it  does  play  a  role  in  shape  and  texture  identification.  Assuming  that  under  static 
loading  the  loading  distribution  is  represented  by  its  Fourier  transform,  fz(u)x,u!y),  then  the 
same  loading  distribution  traveling  along  the  surface  with  a  constant  x-velocity,  vt  will  result 
in  a  loading  distribution  which  has  a  Fourier  transform  /.(wj,  toy)  exp(— juxvtt).  (This  result 
is  easily  obtained  directly  from  the  definition  of  the  Fourier  transform.)  This  means  that  a 
a;,. -comp  one  lit,  of  the  loading  distribution  will  give  rise  to  a  temporal  variations  of  frequency 

ujt  =  vxvt.  (5-1) 

Details  on  this  phenomenon  are  to  be  found  in  Fig.  5-3,  where  we  superimpose  a  typical  spatial 
frequency  distribution  onto  the  observed  temporal  sensitivity  properties  of  human  receptors. 
The  relation  in  Eq.  (5.1)  appears  as  a  line  on  a  loglog  (ut,ux)  plot  and  the  line  is  translated 
sideways  by  changing  the  stroking  velocity  vt.  A  consequence  of  this  result  is  that  the  stroking 
velocity  can  be  selected  based  on  a  desired  spatial  frequency  range  to  match  the  temporal 
sensitivity  range  of  the  sensor/receptor. 


100 


Sensor  Codes 


Sensor  coding  refers  to  how  the  sensors  signals  are  pre-processed  or  what  part  of  its  content 
is  used  for  object  property  identification.  Four  principal  codes,  Cj  —  C4,  have  been  identified 
[34,  40]  based  on  the  assumption  that  we  have  a  distributed  sensors  response  signals  ri(x,y,t), 
that  are  the  response  to  some  stimuli,  where  i  represents  different  sensor  channels. 


1.  Intensive  code  is  when  C  —  Cj  =const.,  i.e.  all  spatial  and  all  temporal  variations  has 
been  integrated  or  averaged  out. 

2.  Temporal  code  is  when  C  =  ^(t),  i.e.  all  spatial  variations  have  been  removed. 

3.  Spatial  code  is  when  C  —  C2(x,y),  i.e.  all  temporal  variations  have  been  removed. 

4.  Spatio-temporal  code  is  when  C  =  C4(x,y,t),  i.e.  at  least  one  channel  preserves  infor¬ 
mation  about  both  spatial  and  temporal  variations. 


All  sensor  channels  may  not  be  used  to  identify  every  object  property  and  it  is  mainly  for 
complicated  object  properties  such  as  texture  that  we  can  expect  more  than  one  channel  to  be 
used.  We  will  now  discuss  how  we  match  the  sensor  channels  and  codes  with  various  object 
properties,  which  we  will  also  define  along  the  way. 


2.1  Identification  of  Object  Properties 


In  this  section  we  will  give  details  on  how  one  can  go  about  identifying  object  properties.  We 
will  go  through  each  of  the  properties  listed  in  Table  5.1,  define  it  and  give  a  modeling  method 
suitable  for  identification,  and  further  delineate  the  appropriate  exploratory  procedures  and 
relevant  sensor  channels  and  codes. 


Scale  Segmentation 


Most  of  object  properties  that  can  be  identified  using  touch  are  a  part  of  the  geometry  of 
the  object  and  one  can  define  different  object  properties  by  segmenting  the  spatial  scale  and 
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Table  5.3:  Lists  definitions  of  geometric  object  properties  and  how  the  different  scale  segments  arise 
from  human  sensing  abilities. 


Name 

Scale 

Comment 

Shape 
Texture 
fi Texture 

Lower  limit  set  by  sensor  resolution 

Lower  limit  set  by  estimated  position  accuracy 
Practical  lower  limit  set  by  sensor  sensitivity 

identifying  a  scale  interval  with  an  object  property.  This  is  necessary  because  there  will  always 
be  an  upper  and  a  lower  limit  on  the  size  of  the  geometric  features  of  interest,  but  real  objects 
will  continue  to  have  geometric  variations  right  down  to  atomic  scale.  We  define  here  three 
object  properties  by  segmenting  the  spatial  scale  and  Figure  5-4  shows  this  schematically,  but 
Table  5.2  in  numbers.  We  name  the  scales  Shape,  Texture  and  p.  Texture,  and  where  the  limits 
are  decided  approximately  by  human  sensing  abilities. 

Having  defined  shape  using  scale  segmentation  we  can  now  simply  define  shape  identification 
as  follows: 


Definition  1  Shape  identification:  Identification  of  the  geometric  boundary  of  a.n  object  to  a 
given  resolution. 


Shape  Identification 


Humans  recognize  common  objects  reliable  in  a  matter  of  few  seconds  [39,  43]  and  based  on 
observation  made  in  human  experiments  the  identification  procedure  has  been  separated  into: 


1.  Scanning,  i.e.  finding  the  object. 

2.  Positioning,  establishing  extents  or  a  point  of  reckoning. 

3.  Enclosure  with  one  or  both  hands  to  get  an  initial  size  estimate. 

4.  Contour  following  where  the  outlines  and  contours  of  the  object  are  traced. 


The  trend  is  therefore  to  begin  the  identification  on  a  coarse  scale  and  the  proceed  to  the  finer 
scales.  This  approach  appears  intuitively  feasible  and  it  is  natural  to  assume  that  the  same 
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approach  could  work  for  robots  as  well.  This  approach  may  however,  not  always  be  optimal 
e.g.  if  the  task  at  hand  is  to  classify  objects  that  only  differ  in  some  detail  at  a  fine  scale, 
then  the  coarse  scale  identification  would  only  give  limited  information.  In  order  to  compare 
or  recognize  an  object  we  must  somehow  be  able  to  register  its  features  and  for  that  we  use 
object  models.  The  problem  we  face  is  to  construct  a  geometrical  model  of  an  object  from 
surface  measurements.  A  similar  problem  occurs  in  machine  vision  and  there,  fortunately  for 
us,  considerable  research  has  been  done  and  several  approaches  to  this  problem  have  been 
developed  (for  review  cf.  [1,  31,  3]).  There  this  problem  has  been  divided  into  five  subtasks: 


1.  Data  acquisition. 

2.  Low-level  processing. 

3.  Object  representation 

4.  Model  Construction. 

5.  Model  matching. 


and  we  note  that  this  subdivision  is  directly  applicable  to  our  problem.  We  have  already  in 
previous  chapters  discussed  subtasks  1  and  2  and  we  go  therefore  directly  to  number  3. 

Object  representation:  Object  representation  refers  to  method  we  use  to  describe  the  shape 
of  the  object.  Many  representations  have  been  developed,  such  as  vertex-  or  edge  based,  surface 
based,  moment  based,  or  volume  based.  We  believe  surface  representation  is  natural,  since  the 
sensor  data  we  have  available  gives  information  about  a  patch  of  the  surface,  i.e.  the  contact 
area,  but  because  of  the  compliant  fingerpad  will  not  give  direct  information  about  either 
vertices  or  edges. 

Model  Construction:  As  an  example  we  will  give  a  short  description  of  a  surface  based 
object  representation  that  can  be  used  in  shape  identification  [50,  78].  It  is  a  physics-based 
approach  in  the  sense  that  it  utilizes  a  3D  elastically  deformable  finite  element  balloon  model 
with  both  membrane  and  plate  properties.  The  spatial  data  then  acts  as  external  forces  pulling 
the  model  surface  towards  where  it  has  been  measured  and  the  internal  membrane/plane  forces 
strive  to  regularize  the  surface  to  be  smooth  and  continuous  in  areas  where  the  surface  shape 
is  unknown. 

Model  Matching:  Once  we  have  a  shape  model  of  the  object  we  can  attempt  to  match  it  with 
a  known  object.  Many  schemes  designed  for  this  task  have  been  developed  in  the  related  fields 
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of  machine  vision  and  artificial  intelligence.  Examples  of  those  are  feature  matching,  where  a 
distance  between  feature  vectors  for  the  shape  model  and  object  is  used  as  a  criteria  to  decide 
which  object  the  shape  model  most  likely  represents  [27],  Other  more  sophisticated  methods 
are  constantly  being  developed  and  we  believe  detailed  discussion  of  those  is  beyond  the  scope 
of  this  thesis. 


As  a  matter  of  completeness  we  list  in  Table  5.4  the  relevant  exploratory  procedures  as  well 
as  the  relevant  sensor  channels  and  codes. 

Table  5.4:  Lists  exploratory  procedures,  sensor  channels  and  sensor  codes  used  for  shape  identification. 


Exploratory 

Sensor 

Sensor 

Procedure 

Channel 

Code 

Static  Touch 

SA1,  RA1 

Spatial 

Stroking 

Strain  gauge  arrays 

Spatio-temporal 

Texture  Identification 


Definition  2  Texture  identification:  Identification  of  the  geometric  features  of  an  object  that 
are  of  finer  scale  than  sensor  resolution  but  that  can  in  principle  be  reconstructed  by  stroking. 


Texture  Measures  As  defined,  texture  could  be  identified  the  same  way  as  shape  (both  are 
in  fact  geometry  of  the  object)  but  we  propose  a  different  identification  method  solely  based 
on  human  experience  and  abilities.  Since  knowing  the  exact  shape  at  such  a  small  scale  is 
of  limited  interest  to  humans  it  is  natural  to  look  for  a  collective  measures  characterizing  the 
shape  variation.  We  will  now  discuss  different  approaches  to  texture  identification  and  their 
characteristics. 

In  the  field  of  neuro-physiology  considerable  research  has  gone  into  identifying  how  texture 
is  coded  in  the  peripheral  nerve  signals  [76,  77,  48,  12,  33,  29].  There  definitions  of  texture, 
roughness  and  measures  thereof  are  widely  varying.  Usually  specimens  used  in  the  experiments 
are  regular  patterns  such  as  dot  patterns,  that  vary  only  in  their  frequency  or  spacing.  The 
measure  applied  is  usually  single  scale,  smooth  to  rough  or  equivalent.  The  research  suffers  from 
inconsistent  definitions  of  texture  and  non-standard  specimens,  but  still  has  some  interesting 
results.  In  [33]  are  presented  results  which  indicate  that  the  tactual  roughness  perception  for 
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the  set  of  specimens  used,  is  based  on  spatial  variation  in  the  SA1  population  response  (which 
in  fact  is  similar  to  am  defined  in  Eq.  (5.3)  below).  The  perceived  roughness  was  further 
observed  to  be  independent  of  stroking  velocity,  but  to  be  proportional  to  the  load  which  is 
consistent  with  our  definition  of  texture  above.  The  problem  is  however  that  most  natural 
textured  surfaces  are  not  as  simple  as  the  specimens  used  and  can  rarely  be  described  using  a 
single  parameter. 

A  related  problem  is  encountered  in  the  characterization  of  surfaces  in  the  field  of  tribology. 
The  most  common  engineering  measure  is  roughness,  Ra,  which  is  primarily  used  along  with  a 
prescribed  machining  process  to  characterize  a  surface  texture  for  machine  parts.  Engineering 
roughness  is  defined  as  the  mean  absolute  deviation  of  a  height  profile  z(x),  from  its  mean: 


A  plethora  of  other  engineering  measures  and  standards  for  the  measure  of  surface  roughness 
exist  [79]  but  all  suffer  to  some  degree  to  from  the  same  shortcomings  as  the  measures  detailed 
above,  i.e.  that  there  are  surfaces  are  identical  in  terms  of  roughness  measures  but  clearly  have 
different  surface  textures. 

The  third  type  of  measures  are  the  scientific  measures  and  they  are  based  on  a  statistical 
approach,  i.e.  the  surface  is  modeled  as  a  two-dimensional  random  process  and  analyzed  as 
such  cf.  [47,  46,  58,  25].  If  the  surface  shape  is  decided  by  the  effects  of  many  random  and 
independent  processes,  we  can  apply  the  Central  Limit  Theorem  which  gives  us  that  the  result 
will  be  a  Gaussian  surface.  Analysis  of  the  joint  Gaussian  probability  distribution  for  the  height, 
the  first  and  the  second  spatial  derivatives,  p(z,z',z")  then  gives  that  such  a  surface  can  be 
characterized  by  only  three  parameters: 


1 .  Root  mean  square  height  of  surface, 


a  =  E  [(z  -  z0y 


1/2 


0.8i?2 


i.e.  the  standard  deviation  of  surface  height. 

2.  Root  mean  square  slope  of  surface, 

crm  =  E  [(z'  -  Zq)2]  7 


(5.2) 


(5.3) 
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3.  Root  mean  square  curvature  of  z, 


Ok 


(5.4) 


From  these  parameters  others  can  be  derived,  such  as  peak  density 

Ds  =  (Ofc/oVn) 


(5.5) 


and  average  wavelength 


A  =  2TUj/<Jm 


(5.6) 


and  many  more.  If  the  surface  is  close  to  being  Gaussian  these  statistical  parameters  can  be  used 
to  identify  and  compare  textures,  e.g.  by  constructing  a  feature  vector  [27].  Not  all  surfaces  are 
Gaussian  and  there  certainly  will  exist  surfaces  with  the  same  (a,  am,  a /.)  but  that  still  cannot 
be  considered  to  have  similar  surface  textures.  This  can  never  be  completely  avoided  as  we 
cannot  completely  describe  an  arbitrary  surface  geometry  using  a  finite  number  of  parameters, 
hence  the  question  is  how  many  parameters  and  which  do  we  need  to  characterize  a  given 
class  of  surfaces.  A  related  issue  is  that  there  will  always  be  some  part  of  the  geometry  that 
cannot  be  sensed  using  tactile  sensors,  namely  sharp  and  deep  minima,  which  give  arise  to  some 
asymmetry.  Here  we  will  not  elaborate  further  on  these  question  but  end  with  the  conclusion 
that  in  order  to  identify  texture  of  surfaces  one  needs  in  general  more  than  one  parameter, 
three  statistical  parameters  for  a  Gaussian  random  surface  but  more  for  a  general  surface. 


Table  5.5:  Lists  Exploratory  procedures,  sensor  channels  and  sensor  codes  used  for  texture  identification. 


Exploratory 

Sensor 

Sensor 

Procedure 

Channel 

Code 

Stroking 

SA1,  RA1 

Strain  gauge  arrays 

Temporal 

Spatio-temporal 

//Texture  Identification 


Definition  3  //  Texture  identification:  Identification  of  geometric  features  of  an  object  that  are 
of  a  finer  scale  than  stroking  resolution  but  can  still  be  detected  as  temporal  variations. 
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Microtexture  differs  from  texture  only  by  scale,  but  the  scale  boundary  between  them  is  set 
by  the  estimated  sensor  properties  of  the  human  hand.  The  sensor  signals  will  be  the  results 
of  the  interaction  of  geometric  variations  that  are  smaller  than  spatial  resolution  and  all  the 
information  can  therefore  be  contained  in  a  temporal  code.  Further,  the  signals  will  typically 
be  high  frequency  so  suitable  sensors  channels  would  be  PC  for  humans  and  accelerometers 
for  robots.  Stroking  velocity  will  effect  the  temporal  frequency  content  of  that  temporal  code, 
higher  velocity  will  translate  the  signal  into  higher  temporal  frequency  range  as  illustrated  in 
Fig.  5-3.  This  in  turn  has  the  consequence  that  if  the  sensors  have  some  optimal  sensitivity 
in  some  specific  frequency  range,  one  can  translate  the  dominant  spatial  frequency,  if  any,  into 
that  range.  For  that  purpose  Eq.  (5.1)  can  be  rewritten  as  v =  A/*  where  A  is  the  average 
wavelength  of  the  surface  defined  in  Eq.  (5.6),  /*  is  the  optimal  temporal  frequency  of  the 
sensor,  and  v*  is  the  corresponding  optimal  stroking  velocity. 

Table  5.6:  Lists  Exploratory  procedures,  sensor  channels  and  sensor  codes  used  for  //.texture  identifica¬ 
tion. 


Exploratory 

Procedure 

Sensor 

Channel 

Sensor 

Code 

Stroking 

PC 

Accelerometers 

Intensive 

Temporal 

Compliance  Identification 

Definition  4  Compliance  Identification:  the  estimation  of  the  effective  compliance  of  an  object, 
at  the  contact  location,  averaged  over  the  contact  area,  assuming  that  the  object  is  smooth  and 
made  of  homogeneous  and  linearly  elastic  material. 


The  compliance  of  an  object  is  a  material  property  that  refers  to  how  much  the  material 
deforms  with  increased  loading.  The  identified  compliance  should  ideally  not  depend  on  the 
loading  or  geometry  of  the  object.  That  is  however  not  practical  as  it  is  generally  not  possible 
to  distinguish  between  the  effects  of  geometry  and  material  properties  if  both  are  unknown. 
Here  we  will  therefore  neglect  the  effects  of  geometry  and  assume  that  the  object  is  smooth  over 
the  contact  area.  We  are  then  calculating  the  effective  compliance  of  that  surface  averaged  over 
the  contact  area.  When  the  object  is  smooth  and  slowly  varying,  the  contact  area  will  have 
the  shape  of  an  ellipse  and  the  loading  distribution  will  be  ellipsoidal,  i.e.  Hertzian  contact 
conditions  [36].  Further,  when  a  compliant  material  is  contacted  with  compliant  fingerpad 
the  compliances  of  add  the  same  way  as  springs  in  series.  Assuming  that  the  compliance 
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of  the  fingerpad  is  known  we  can  then  infer  the  compliance  of  an  object  from  the  combined 
compliance.  For  a  linearly  elastic  material  the  compliance  depends  on  two  constants,  E ,  v, 
the  Young’s  modulus  and  Poisson’s  ratio.  How  they  combine  depends  on  what  exploration 
procedure  is  used,  indentation  or  stretch,  or  equivalently  the  direction  of  the  load,  normal  or 
tangential,  but  the  general  form  is  similar. 


Normal  Compliance:  In  [36]  it  is  shown  that  if  the  contact  area  has  the  form  of  the  circle 
with  radius,  a,  then  the  overall  normal  displacement,  8Z ,  is  related  to  the  total  normal  load,  Fz 

by 

6.  =  — — .  (5.7) 

*  T7I*  A  ..  v  ' 


E*  4  a 


Where  1/E*  is  the  combined  normal  compliance,  given  by 


1  1-^  l~vl 

E*  Ef  E0 


(5.8) 


and  where  the  subscript,  /  refers  to  the  material  properties  of  the  finger,  but  o  to  that  of  the 
object.  This  relation  is  also  valid  if  linear  viscous  effects  are  present,  [30,  44,  36],  but  then 
only  while  the  contact  area  is  increasing2  a  >  0,  but  then  E*  =  E*(s)  is  in  the  form  of  a 
Laplace  transfer  function  describing  the  rate  dependent  behavior  of  the  material.  In  this  case 
one  can  simply  apply  system  identification  techniques  from  control  theory  (cf.  [45])  to  identify 
the  parameters  of  E*(s). 


In  the  more  general  case,  when  the  contact  area  is  an  ellipse,  with  axes  a  and  b,  the  relation 
in  Eq.  (5.7)  is  only  slightly  modified  and  becomes 


8>z  — 


1  3F:  bV(  \ 

E*  4(ab)V*a  {G> 


(5.9) 


where  A'(e)  is  a  complete  elliptic  integral,  and  e  =  (1  —  (6/a)2)1/2,  This  solution  is  similarly 
valid  for  viscoelastic  materials  while  the  contact  area  is  non-decreasing  [85]. 


Compliance  by  Tapping:  It  has  been  observed  that  humans  can  quite  effectively  discriminate 
between  differently  compliant  materials  by  tapping  it  with  a  stylus.  Sound  is  evidentially  an 
important  cue  and  besides  that  experiments  also  indicate  that  the  initial  rise  time  of  the  loading 
transmitted  to  the  fingerpad  is  the  most  important  tactile  feature  used  [41]. 

“This  solution  becomes  invalid  if  d(f)  at  any  time  becomes  negative  and  then  the  solution  becomes  significantly 
more  complicated,  which  form  depends  on  the  loading  history  [82,  8]. 
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Table  5.7:  Lists  Exploratory  procedures,  sensor  channels  and  sensor  codes  used  for  normal  compliance 
identification. 


Exploratory 

Sensor 

Sensor 

Procedure 

Channel 

Code 

Indenting 

SA1,  SA2 

Intensive 

Strain  gauge  array 

Spatial 

Force  sensors 

Spatio-temporal 

These  results  can  be  explained  theoretically  by  calculating  the  time  to  peak  force,  f*,  for  a 
stylus  that  impacts  a  linearly  elastic  half-space  [49,  36].  The  time  depends  on  the  shape  of  the 
end  of  the  stylus  but  for  spherical  and  flat  ends  t*  can  be  calculated  to  be 


t 


Spherical 


l-I/2 

E 


2/5 


+*  _ 

rFlat  — 


l-!/2 


1  1/2 


E 


(5.10) 


Figure  5-5  further  shows  the  response  force  trajectory  for  three  different  materials,  from  which 
we  conclude  that  this  relation  between  the  time  to  peak  force  and  compliance  sufficiently  ex¬ 
plains  experimental  results. 


Table  5.8:  Lists  Exploratory  procedures,  sensor  channels  and  sensor  codes  used  for  normal  compliance 
identification  using  tapping. 


Exploratory 

Procedure 

Sensor 

Channel 

Sensor 

Code 

Tapping 

RAl,  PC 

Force  sensors 
Accelerometers 

Temporal 

Tangential  Compliance:  The  lateral  behavior  is  slightly  different  because  there  the  effects 
of  friction  must  be  included.  When  simultaneous  normal  and  tangential  loads  are  applied 
then  there  is  a  central  region,  homothetic  with  the  contact  area,  where  no  slip  occurs,  and  a 
surrounding  region  with  microslip  [36,  14]  (see  also  section  4,  chapter  3).  When  the  contact 
region  is  circular,  with  radius  a,  the  relation  between  the  tangential  displacement,  Sx,  and  the 
tangential  load,  Fx  is 

c  _ 

x  Ef  4  a 


-(-A) 


2/3 


(5.11) 


Where  1  }E\  is  the  combined  lateral  compliance,  given  by 


J_  _  (2  -  vf)(l  +  uf)  (2  -  v0)(l  +  v0) 
Ef  2  Ef  2  E0 
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where  the  subscript,  /  refers  to  the  material  properties  of  the  finger,  but  o  to  that  of  the  object. 
We  note  that  if  Ef  =  Ea  and  Uf  =  v0,  then  the  ratio  of  normal  compliance  to  lateral  compliance 
is  2(1  —  v)/(2  —  v)  which  goes  from  1  to  1.5  as  u  goes  from  0  to  0.5,  i.e.  lateral  and  normal 
compliances  are  comparable  in  this  case. 


For  non-axisymmetric  objects,  i.e.  when  the  contact  area  is  an  ellipse,  similar  relations  have 
been  derived  and  similarly  there  is  a  correction  factor  a  function  of  the  eccentricity,  e,  but  the 
relations  are  not  as  simple  as  above  and  we  confer  the  reader  to  [14]  for  details. 


Table  5.9:  Lists  Exploratory  procedures,  sensor  channels  and  sensor  codes  used  for  lateral  compliance 
identification. 


Exploratory 

Sensor 

Sensor 

Procedure 

Channel 

Code 

Stretch 

SA1,  SA2 

Intensive 

Static  touch 

Strain  gauge  array 

Spatial 

Force  sensors 

Spatio-temporal 

Layered  Compliance:  A  variant  of  compliance  identification  is  the  case  when  rve  have  an 
object  that  is  composed  of  different  layers  that  can  slide  relative  to  each  other,  e.g.  skin 
over  a  bone  or  a  muscle.  The  task  is  then  to  detect  the  shape  and/or  compliance  of  the 
deeper  layer,  e.g.  finding  lumps  in  tissue  which  is  an  important  application  of  tactile  sensing  in 
medical  practice  [48].  One  way  of  achieving  that  is  to  modify  the  transfer  functions  used  in  the 
decoding  algorithms  (see  chapters  3  and  4)  by  increasing  the  parameter  for  sensor  depth  so  that 
the  signals  estimates,  the  output,  are  the  estimated  displacements  and  loading  distributions  at 
the  depth  of  the  second  layer.  Its  shape  can  then  be  identified  and  its  compliance  can  also  be 
estimated  using  the  same  relations  as  above. 

Table  5.10:  Lists  Exploratory  procedures,  sensor  channels  and  sensor  codes  used  for  layered  compliance 
identification. 


Exploratory 

Sensor 

Sensor 

Procedure 

Channel 

Code 

Stretch 

SA1,  SA2 

Intensive 

Stroking 

Strain  gauge  array 

Spatial 

Force  sensors 

Spatio-temporal 
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2.2  Identification  of  the  State  of  Contact 


Definition  5  State- of- Contact  Identification:  the  identification  of  the  nature  or  mode  of  the 
interaction  between  the  object  and  the  fingerpad  from  sensor  signals  and  estimation  of  the  related 
parameters. 


As  stated  in  the  definition,  then  the  concept  state-of-contact  refers  to  which  mode  the 
interaction  between  the  object  and  the  fingerpad  is  at  the  contact  area.  Is  there  adhesion, 
sliding  with  friction  or  something  in-between?  State  of  contact  differs  from  the  object  properties 
delineated  above,  in  that  it  is  not  strictly  a  property  of  the  object  only,  but  is  rather  an  interface 
property.  For  example  friction,  which  is  very  important  in  tactile  sensing  depends  in  reality 
on  the  load,  sliding  velocity,  temperature,  and  surface  texture  in  addition  to  individual  and 
combined  material  properties  [56,  4].  Identifying  friction  parameters,  e.g.  the  friction  coefficient 
Up,  therefore  has  little  meaning  except  as  an  interface  property  at  that  location  and  instant 
only. 

Adhesion:  Adhesion  is  mainly  important  for  very  soft  and  smooth  objects3  such  as  gels  or 
gel-like  materials.  A  clear  cue  for  adhesion  is  obviously  negative  pressure  in  the  contact  area.  It 
is  however  unknown  if  adhesion  plays  a  role  in  human  haptics  as  fingerpad  irregularities,  such 
as  fingerprints  prevent  in  most  cases  adhesion.  Instead  the  finger  ridges  give  better  hold  in  wet 
and  lubricated  contact  conditions,  the  same  principle  as  is  behind  tiretracks. 

Microslip:  In  section  4.  chapter  3  we  discussed  in  detail  what  happens  when  tangential  loading 
is  added  to  an  object  that  is  already  under  normal  loading,  assuming  that  Coulomb’s  law  and 
axisymmetric  Hertzian  conditions  are  valid.  For  sake  of  completeness  we  will  review  those 
results  here  but  refer  the  reader  to  chapter  3  or  [36]  for  further  details. 

As  one  increases  the  tangential  load,  Ft  from  zero  it  turns  out  that  slip  does  not  occur 
instantaneously  when  the  friction  limit,  Ft  —  pjFz  is  reached.  The  object  sticks  to  the  finger¬ 
pad  in  a  circular  central  region  with  radius,  c,  but  microslip  occurs  in  the  surrounding  region 
starting  from  the  edge  of  the  contact.  In  the  microslip  region  the  tangential  loading  is  related  to 
the  normal  loading  through  Coulomb’s  law  but  in  the  sticking  region  it  follows  load  distribution 
giving  uniform  displacement.  Figure  5-6  illustrates  this  phenomena,  showing  the  loading  situ¬ 
ation  in  Fig.  5-6a,  the  stick  and  microslip  areas  superimposed  on  a  grid  showing  the  resulting- 
tangential  displacement  in  Fig.  5-6c.  Profiles  of  the  tangential  and  normal  load  distributions 

3We  exclude  sticky  films,  glues  etc.  from  our  discussion  here. 
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are  shown  in  Fig.  5-6c  and  finally  Fig.  5-6d  shows  the  relation  between  the  applied  tangential 
load  and  the  stick  radius,  which  turns  out  to  be  given  by 

-%-  =  !-  (c/af.  (5.13) 

Vfi'z 

Since  the  tangential  loading  distribution  has  a  sharp  edge  at  the  outer  limits  of  the  stick  region 
(r  =  c)  its  extents  can  be  inferred  from  the  tangential  load  distribution  and  along  with  the 
contact  radius  used  to  calculate  the  friction  coefficient  from  Eq.  (5.13)  without  any  global  slip 
of  the  contact.  This  is  a  very  valuable  information  that  can  be  used  to  predict  the  onset  of 
slip  before  it  happens.  More  complicated  but  qualitatively  similar  behavior  can  be  expected 
for  more  generally  shaped  object  but  will  not  be  discussed  further  here,  (cf.  [67]). 

Sliding:  If  the  object  is  smooth  so  that  Coulomb’s  law  and  Hertzian  conditions  are  valid  then 
during  sliding  the  tangential  and  normal  load  distribution  will  be  proportional  to  each  other.  In 
the  absence  of  other  spatial  information,  such  as  texture,  the  sliding  velocity  cannot  be  inferred 
from  the  sensors  as  /Texture  will  be  based  on  a  temporal  code  which  cannot  uniquely  decide  the 
sliding  velocity.  Hence,  either  kinesthetic  information,  texture,  or  shape  variations  that  have 
smaller  spatial  wavelength  than  the  fingerpad  are  needed  to  infer  the  sliding  velocity. 


3.  Finger  Control 


In  this  section  we  will  discuss  what  happens  inside  the  Finger  Control  block  in  Fig.  5-1  and  in 
Fig.  5-7  we  have  expanded  that  very  block  and  shown  its  inside  structure.  In  Fig.  5-7  we  have 
divided  the  controller  into  two  main  paths,  tangential  controller  and  normal  controller.  We  will 
first  discuss  the  tangential  controller  as  it  is  more  involved  and  mainly  consider  stroking.  Other 
exploratory  procedures  are,  from  the  viewpoint  of  control,  much  simpler  and  can  as  such  be 
considered  special  cases  of  stroking. 


3.1  Tangential  Control 


During  stroking  the  tangential  controller  must  ensure  that  the  fingerpad  maintains  contact  and 
follows  the  shape  of  the  object,  while  maintaining  the  contact  area  within  the  extents  of  the 
fingerpad.  As  discussed  in  section  1.  we  assume  here  that  the  desired  direction  of  stroking,  i.e. 
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the  tangential  velocity  vector,  is  provided  by  an  outside  source  (supervisor).  However,  in  order 
to  keep  the  contact  region  within  the  fingerpad  area  during  stroking,  the  fingerpad  will  need 
to  rotate  as  it  moves  over  the  surface.  How  this  simultaneous  rolling  and  sliding  is  combined  is 
decided  by  the  contact  or  stroking  kinematics. 


Stroking  Kinematics 


Equations  describing  the  kinematics  of  an  object  and  a  finger  in  contact,  subjected  to  simulta¬ 
neous  relative  rolling  and  sliding  have  been  available  in  the  literature  for  some  time  [53].  They 
have  recently  been  extended  to  include  multiple  fingers  and  contact  dynamics  [57,  55].  Here  we 
will  focus  on  the  special  and  relatively  simple  case  of  stroking,  i.e.  when  the  center  of  contact 
does  not  move  relative  to  the  fingerpad.  Figure  5-8  shows  a  fingerpad  stroking  an  object  and 
defines  a  coordinate  system  fixed  in  the  fingerpad. 

Contact  kinematics  then  demand  that  the  translation  velocities  vx,  vy )  and  the  rotational 
velocities,  ( 6X ,  0y )  of  the  fingerpad  be  related  through  the  curvature  KQb}  of  the  object  at  the 
contact  location  as  follows 


Vx 

II 

1 - " 

H 

.  ~vy  . 

When  the  fingerpad  is  flat4  A'obj  is  a  2x2  matrix  of  second  spatial  derivatives  of  the  objects 
shape  taken  with  respect  to  the  coordinated  system  defined  in  Fig.  5-8.  The  relation  in  Eq. 
(5.14)  is  obtained  under  the  assumption  that  the  object’s  shape,  as  well  as  its  first  and  second 
spatial  derivatives  vary  smoothly  in  the  contact  area.  Even  if  that  is  true,  measurements  of 
derivatives  are  always  tricky,  especially  in  the  presence  of  imperfect  measurements.  Addition¬ 
ally,  for  our  purpose,  i.e.  stroking,  we  are  only  interested  in  the  average  or  overall  curvature  of 
the  contact  area. 


Overall  Curvature  Estimation 


The  problem  of  finding  the  overall  curvature  of  the  contact  area  can  equivalently  be  stated  as 
fitting  a  constant  curvature  displacement  distribution  to  the  estimated  surface  displacements. 
The  curvature  should  only  be  estimated  within  the  contact  area,  defined  here  as  being  where 

4Similar  relations  hold  when  the  fingerpad  is  also  curved,  cf.  [53]  for  details. 
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the  normal  load  is  positive  (i.e.  indicating  pressure)  and  therefore  we  also  need  the  estimated 
normal  load  distribution.  An  estimated  distribution  can  be  fitted  to  a  constant  curvature 
distribution  using  standard  fitting  techniques.  Fitting  in  the  least  square  sense  can  for  example 
be  obtained  by  matching  the  moments  of  the  distributions  and  use  those  relations  to  solve 
for  the  parameters  of  the  constant  curvature  distribution.  A  constant  curvature  displacement 
distribution  can  generally  be  described  using  six  parameters 

uZo  —6  —  +  yk'2  CXU  +  Dx  +  Ey'j  (5.15) 

We  observe  that  the  zeroth,  first  and  second  moments  will  give  all  six  parameters,  but  the 
moments  are  defined  as: 


/o  =  fcuzodxdy 

lx  =  fc  xuzodxdy  Iy  =  fc  yuzodxdy  (5,16) 

I xx  =  fc  x2uzodxdy  Ixy  =  fc  xyuZodxdy  Iyy  =  fc  y2uzodxdy 

All  the  integrals  are  taken  over  the  contact  area  C  which  is  known  from  the  surface  load 
distribution  as  discussed  above. 


The  entries  in  Ka bj  are  then  finally  obtained  by  calculating  each  of  the  moments  for  both 
the  constant  curvature  distribution  in  Eq.  (5.15)  and  the  estimated  displacement  distribution 
U‘s  and  then  solving  for  for  the  .4,  B  and  C  values  which  are  the  values  needed,  since 


Aobj 


A  C 
C  B 


When  viscous  effects  are  negligible  and  quasi-static  contact  conditions  are  valid,  we  know 
that  a  Hertzian  (i.e.  ellipsoidal)  load  distribution  will  give  a  constant  curvature  displacement 
distribution.  We  can  therefore  work  directly  with  the  load  distribution,  i.e  fit.  the  estimated 
normal  load  distribution  to  a  Hertzian  distribution,  and  hence  avoid  the  problem  of  explicitly 
estimating  the  area  of  contact.  Hertzian  load  distribution  can  generally  also  be  described  using 
six  parameters,  e.g. 

fz0  =  4>z\J 1  -  (Afx2  +  Bfy 2  +  Cfxy  +  Dfx  +  Efy)  (5.17) 

Similarly  as  above,  the  parameters  of  this  relation  can  be  estimated,  in  the  least  square  sense, 
by  matching  the  zeroth,  first  and  second  moments  of  the  estimated  and  Hertzian  distributions. 
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The  entries  for  K0bi  can  then  be  obtained  from  known  Hertzian  load-displacement  relations, 
cf.  [36]. 

Example: 

As  an  example  we  now  show  how  the  method  delineated  above  can  be  used  to  estimate  the 
overall  curvature  of  a  textured  surface.  Figure  5-9a  shows  the  surface  from  a  perspective  view 
but  Fig.  5-9c  from  above  as  a  height-density  plot.  Figure  5-9b  shows  that  part  of  the  surface 
that  contacts  the  fingerpad  at  this  particular  instant  and  similarly  does  Fig.  5-9d  show  the  same 
from  above  as  height-density  plot.  We  note  that  the  contact  area  is  not  a  simply  connected 
region  and  that  one  cannot  reliably  estimate  the  curvature  from  single  location  within  that 
area. 


As  discussed  above  we  use  the  six  moments  defined  in  Eq.  (5.16)  to  obtain  six  equations 
that  we  use  to  solve  for  the  parameters  of  interest  in  the  constant  curvature  defined  in  Eq. 
(5.15).  We  combine  this  relations  in  a  matrix  equation 


where,  for  example, 


and  further 
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(5.18) 


where  as  before  u~0  is  the  estimated  surface  displacement. 


Figure  5-10  shows  the  results  of  the  above  applied  to  the  surface  show  in  Fig.  5-9.  Figure 
5-10a  shows  it  as  a  surface  plot  with  the  fingerpad  surface  displacement  superimposed  onto  the 
resulting  constant  curvature  surface  and  Fig.  5-10b  shows  the  same  but  from  a  profile  view.  We 
note  that  the  constant  curvature  surfaces  captures  the  overall  curvature  of  the  surface  despite 
the  widely  varying  curvature  of  the  texture. 
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Fitting  the  estimated  displacement  distributions  this  way  gives  us  an  easily  calculable  and 
the  estimate  is  robust  since  we  avoid  taking  derivatives.  However,  care  must  be  taken  when 
the  contact  area  becomes  very  small.  The  loads  will  the  become  concentrated  and  the  effects 
of  discrete  sensors  can  then  lead  to  inaccurate  estimates  of  the  overall  curvature.  It  is  however 
a  different  question  how  we  should  react  to  sharp  points  or  corners  and  otherwise  regulate  the 
normal  load  so  that  the  contact  area  is  large  enough  but  does  not  lead  to  overly  high  load 
concentrations  that  can  damage  the  sensors. 


3.2  Normal  Control 


The  problem  of  maintaining  a  suitable  normal  force  while  moving  along  a  surface  has  been 
under  investigation  in  the  field  of  robotics  for  some  time  [2]  where  the  term  compliant  motion 
control  has  been  used.  Solution  to  this  problem  have  been  formulated  in  a  number  of  ways  using 
different  control  schemes.  The  most  straightforward  approach,  in  our  case,  is  to  directly  control 
the  total  normal  force,  i.e.  a  force  control  scheme  in  the  normal  direction.  A  variant  of  this 
approach  would  be  to  control  or  regulate  the  maximum  normal  pressure  or  some  other  norm  or 
measure  of  the  normal  load  distribution.  Other  control  schemes,  such  as  impedance  control  can 
also  be  applied,  but  since  we  already  have  measurements  of  the  contact  force  they  may  not  be 
needed  but  one  of  its  strengths  is  that  force  measurements  are  not  needed  [2,  26].  Impedance 
control  might  however  prove  useful  during  the  phases  when  the  hand  is  in  free  motion  prior  to 
or  during  contact  [11], 


4.  Conclusions 


In  this  chapter  we  have  discussed  the  computational  aspects  of  exploration  of  objects  and 
environments.  We  have  analyzed  exploration  strategies  in  terms  of  the  exploration  procedures, 
sensor  channels,  codes  and  object  properties  that  are  available  using  tactile  sensing.  In  addition 
to  providing  a  general  and  common  structure  positioning  exploration  of  objects  relative  to  other 
aspects  of  tactile  sensing,  we  have  analyzed  each  of  the  above  components  in  detail  and  explicitly 
stated  related  formulas  where  appropriate.  We  have  finally  discussed  the  implementation  of 
proposed  exploration  strategies  in  terms  of  different  control  schemes. 

A  framework  that  organizes  the  different  aspects  that  relate  to  computational  theory  of 
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haptics  is  a  very  important  milestone  -  if  it  gains  acceptance  among  researchers.  We  still  must 
recognize,  that  we  cannot  claim  that  results  presented  here  are  based  on  the  only  viable  or 
reasonable  approach.  For  example,  the  intermediate  step  of  surface  signal  identification,  where 
we  explicitly  estimate  both  the  surface  load  and  surface  displacement  distributions,  may  not 
be  necessary.  Instead,  the  object  properties  could  be  identified  directly  from  the  sensor  signals. 
We  however  believe  that  an  explicit  formulation  like  ours  is  needed  first,  and  that  our  work  has 
value  in  being  a  reasonable  starting  point  and  providing  understanding  useful  for  future  work. 
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Figure  5-1:  Shows  a  diagram  that  explains  the  structure  and  the  relation  between  the  various  problems 
that  constitute  a  computational  theory  of  haptics.  The  Surface  Signal  Identification  problem  was  treated 
in  Chapters  3  and  4  but  in  this  chapter  we  will  discuss  the  Exploration  Strategy,  and  Finger  Controller 
blocks.  The  Exploration  Strategy  block  refers  to  the  problem  of  inferring  object  properties  such  as 
shape,  texture  and  compliance  or  state  of  contact  from  the  estimated  surface  signals.  This  involves 
determining,  given  a  property  that  is  to  be  explored,  what  kind  of  sensor  information  and  what  kind  of 
action  is  needed  for  successful  exploration.  The  Supervisor  shown  in  the  figure  represents  an  external 
source  that  supplies  additional  information  required  to  identify  each  object,  e.g.  the  stroking  direction 
for  shape  identification.  The  Finger  Controller  block  represents  the  calculation  of  tracking  trajectories 
for  the  fingerpad  from  a  given  desired  action,  local  shape,  and  contact  force  distribution  information. 
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Temporal  frequency,  cot 


Figure  5-2:  Schematically  shows  the  temporal  and  spatial  frequency  characteristics  of  human  receptors 
as  ranges  of  sensitivity. 
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co./ 271  [Hz] 


Figure  5-3:  Shows  a  typical  spatial  frequency  distribution  superimposed  onto  the  temporal  sensitivity 
properties  of  human  receptors.  The  relation  in  Eq.  (5.1)  appears  as  a  line  on  the  loglog  (u>t ,  ux  )  plot 
(upper-right)  and  the  line  translates  sideways  by  changing  the  stroking  velocity  vt ■ 


Original  Geometry 


Texture  jlTexture 
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Figure  5-4:  Shows  scale  segmentation  schematically  and  how  the  original  geometry  is  composed  of 
Shape,  Texture  and  //Texture. 
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Impact  Force  for  Different  Stiffnesses 


Figure  5-5:  Impact  response  trajectories  for  differently  compliant  materials.  The  trajectories  are  nor¬ 
malized  with  respect  to  that  of  the  intermediate  compliant  material.  We  observe  that  there  is  a  direct 
relation  between  the  time  to  peak  force  and  material  compliance,  as  indicated  in  Eq.  (5.10),  and  that  it 
can  be  used  for  material  discrimination. 
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(a)  Configuration  (b)  Tangential  Displacement 


Figure  5-6:  Illustrates  the  behavior  of  axisymmetric  Hertzian  solid  under  simultaneous  normal  and 
tangential  load,  (a)  shows  the  loading  configuration,  (b)  shows  the  stick  and  microslip  region  superim¬ 
posed  onto  a  grid  showing  the  resulting  tangential  displacement,  (c)  shows  profiles  of  the  tangential  and 
normal  load  distributions  and  finally  (d)  shows  the  relation  between  normal  load,  tangential  load,  stick 
radius,  contact  radius  and  friction  coefficient,  i.e.  Eq.  (5.13). 
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Figure  5-7:  Shows  a  diagram  that  expands  the  Finger  Controller-block  in  Fig.  5-1.  There  are  two  main 
paths:  i)  tangential  control  and  ii)  normal  control.  The  tangential  control  takes  as  inputs  estimated 
surface  signals  and  stroking  direction  and  velocity  from  a  supervisor.  It  gives  as  output  the  desired 
tangential  translation  and  rotation  velocities  ( vx ,  vy , 8X ,  6y)  after  going  through  the  steps  of  Curvature 
Estimation  and  Stroking  Kinematics.  The  normal  control  path  takes  as  input  the  estimated  surface  load 
distribution  and  gives  as  output  desired  total  normal  force.  Both  normal  and  tangential  outputs  are  fed 
into  a  Low-level  Controller  which  we  assume  is  capable  of  tracking  the  desired  signals  with  sufficient 
accuracy. 
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Figure  5-8:  Shows  a  fingerpad  stroking  an  object  and  defines  a  coordinate  system  fixed  in  the  fingerpad. 
With  the  tangential  velocities  (vx,  vy)  and  rotational  velocities  about  the  respective  axes  (8X,  6y)  . 
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Figure  5-9:  Shows  a  textured  surface  and  the  part  of  the  surface  that  is  in  contact  with  the  fingerpad. 
Figure  (a)  shows  the  actual  surface  from  a  perspective  viewpoint,  (c)  from  above  as  a  height-density 
plot.  Similarly,  shows  figure  (b)  the  part  in  contact  from  a  perspective  but  (d)  from  above.  We  note 
that  the  contact  region  is  not  a  simply  connected  region,  as  there  are  areas  where  there  is  a  gap  between 
the  object  and  the  fingerpad. 
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Figure  5-10:  Shows  a  constant  curvature  surface  superimposed  on  the  texture  surface  in  Fig.  5-9 
obtained  by  matching  the  moments,  defined  in  Eq.  (5.16),  of  the  constant  curvature  surface  and  the 
surface  displacement  within  the  contact  area.  Figure  (a)  shows  it  as  a  surface  plot  and  (b)  shows  it 
from  a  profile  view  where  the  contact  region  is  marked  with  thickened  lines.  We  note  that  the  constant 
curvature  surface  captures  the  overall  curvature  of  the  surface  despite  the  widely  varying  curvature  of 
the  texture. 
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Chapter  6 


Concluding  Remarks  and 
Suggestions  for  Future  Work 


1.  Conclusions 


In  this  thesis  we  have  presented  a  basic  computational  theory  of  haptics.  We  divide  the  problem 
into  Identification  and  Control  problems.  We  further  separate  the  identification  problem  into 
Surface  Signal  Identification  and  Object  Property  Identification  problems  as  shown  in  Fig.  5- 

1.  We  have  analyzed  each  of  these  problems  and  presented  our  results  in  four,  by  and  large, 
independent  chapters,  each  with  their  concluding  remarks,  but  they  can  be  summarized  as 
follows: 


1.  In  order  to  correctly  decode  a  general  3D  tactile  signal  from  a  sensor  array,  we  need  a 
minimum  of  three  independent  sensor  signals,  and  we  show  that  (pn,  exz ,  eyz)  is  the  only 
stress/strain  signal  combination  without  a  directional  bias. 

2.  We  give  solutions  that  are  valid  under  static  or  quasi-static  conditions  as  well  as  ones  that 
include  the  effects  of  viscoelasticity  and  inertia  and  further  give  conditions  under  which 
the  inertia  can  be  ignored. 

3.  The  transfer  function  approach  developed  in  the  thesis  is  shown  to  allow  us  to  infer 
qualitative  properties  and  behavior  of  the  solution  without  specifying  the  inputs.  The 
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transfer  function  approach  further  gives  us  access  to  powerful  signal  processing  tools, 
such  as  the  Wiener  and  Kalman  filter  which  we  evoke  to  solve  the  decoding  problem.  The 
decoding  problem  being  the  problem  of  reconstructing  the  surface  signals  from  sensor 
signals,  in  a  robust  manner,  in  the  presence  of  sensor  noise. 

4.  We  analyze  exploration  strategies  in  terms  of  exploration  procedures,  sensor  channels, 
codes  and  the  object  properties  that  are  available  using  tactile  sensing  and  explicitly 
state  related  formulas  where  appropriate. 

5.  Finally,  we  discuss  implementation  issues  in  terms  of  different  control  schemes. 


The  above  results  are  valid  for  homogeneous  and  linearly  elastic  materials  and  while  in- 
finitismal  displacement  assumptions  hold.  These  are  idealized  conditions  and  it  is  therefore 
important  to  consider  what  happens  when  they  are  not  completely  satisfied.  The  tools  and 
methods  we  apply  have,  however,  been  applied  with  good  results  to  various  problems  where 
those  assumptions  may  not  always  hold,  most  notably  the  methods  developed  around  linear 
elasticity  in  mechanics  and  also  the  Kalman  filter  in  navigation  and  control.  We  however 
cannot  state  unilaterally  that  the  solutions  given  in  this  thesis  will  not  fail  under  non-ideal 
conditions,  but  that  can  most  likely  only  be  assessed  on  case-by-case  basis  and  then  with  the 
aid  of  experiments. 


2.  Suggestions  for  Future  Work 


The  results  presented  here  essentially  provide  a  skeleton  of  a  computational  theory  or  an  outline 
to  which  details  need  to  be  added.  The  encoding  solutions  are  given  for  the  infinite  halfspace 
but  solutions  for  other  geometries,  such  as  layered,  finite  thickness,  finite  size,  cylindrical  or 
spherical,  might  more  accurately  model  real  fingerpads  and  therefore  prove  more  useful. 

The  decoding  solution  given  here  assumes  that  the  sensor  signal  is  linearly  related  to  the 
components  of  the  stress  or  strain  tensor,  but  we  might  want  to  consider  other  nonlinear  mea¬ 
sures  such  as  strain  energy  density  or  maximum  compressive  strain  that  experimental  results 
indicate  might  constitute  the  relevant  sensor  stimulus  in  the  human  tactile  system  [13,  75].  In 
appendix  F  we  state  the  decoding  problem  in  general  terms  where  we  observe  that  the  problem 
is  then  expressed  in  conditions  regarding  the  integrability  of  the  relevant  PDEs.  The  integra- 
bility  of  general  PDEs  is  an  advanced  mathematical  problem  beyond  the  scope  of  this  thesis 
and  we  do  not  discuss  that  further  here  but  refer  the  reader  to  [64]. 
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The  object  identification  problem  has  many  possibilities  for  further  work,  such  as  experi¬ 
mental  studies  of  human  subject  where  the  experimental  paradigms  and  results  are  stated  in 
terms  of  the  concepts  we  have  developed.  As  an  example  one  could  test  the  perception  of 
texture  in  the  terms  of  roughness  measures  or  test  the  hypothesis  that  humans  perceive  texture 
and  //texture  differently,  the  former  being  independent  of  stroking  velocity. 
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Appendix  A 


Derivation  of  Static  Transfer 
Function  Matrices 


We  start  with  the  differential  equations  of  equilibrium1 

=  0  hj  -  +  (A-!) 

and  the  linear  stress  strain  relations  according  to  Hooke’s  law 

(T ij  —  \Sij6ij  +  2fl£ij  i,j  —  (A-2) 


where 


£ij  =  ~  (Uij  +  Uj'i )  i,j  =  x,  y ,  # 


and  Sij  =  1  if  and  only  if  i  =  j,  and  is  zero  otherwise.  Further,  the  Lame  constants  A  and  p  (shear 
modulus)  are  in  terms  of  Young’s  modulus  and  Poisson  ratio: 


vE  E 

(l  +  i/)(l  —  2v) '  2(l  +  i/)‘ 


'Here  we  use  tensor  notation  for  compactness,  where  a  ”,  j"  in  subscript  means  a  partial  derivative  with 
respect  to  the  independent  variable  represented  by  j  and  repeated  symbols  in  subscript  stand  for  summation 
over  all  indices. 
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In  order  to  solve  these  nine  partial  differential  equations  we  make  use  of  the  two  dimensional  Fourier 
transform  defined  by  [5] 


r  r  +  OO 

f(ux,uy)  =  //  f(x,y)e~J^:X+u,'-'y'>dxdy 

J  J  —  OO 


i  r  r +°° 

f(x,y)  =  dwxdo;B. 


Assuming  that  all  integrals  are  well  behaved  as  | a;,  y\  — >  oo,  we  apply  this  transform  to  each  term  in  the 
equations  in  (A.l)  and  (A. 2).  Then  the  operator  is  replaced  by  multiplication  by  jux .  and  #  by 
juix.  We  further  use  D  to  notate  The  9  partial  differential  equations  have  now  been  transformed 
into  ordinary  differential  equations  in  D  ( becomes  ^).  We  write  (A.l)  and  (A. 2)  using  matrix-vector 
notation  as 


and 


K Ta  =  0,  7  =  Ku, 


cr  —  C  7 


where 


]ux 

0 

0 


0 

JOJy 

0 


0  jU)y 

0  juix 
D  0 


D 

0 


0 

D 


3  & x  3  My  J 


c  = 


A  -f-  2  n 
A 
A 
0 
0 
0 


A 

A  +  2 
A 
0 
0 
0 


A  0 

A  0 

A  T  2  {i  0 

0  p 

0  0 

0  0 


0 

0 

0 

0 

0 

II 


<?xz,  &yz]T ,  and  7  =  [exx,  Eyy,  ezz,  7 xy,  yxz,  yyz]T2.  Eliminating  a  and  E  we  get 


KtCI<u  =  0 


(A.3) 


where  u  =  \ux,  uy,  uz]  and 


ktci<  = 


D2  —  Q,2  —  oj2(02  —  1) 

-U)xLOy{02  -  1) 

j^x{/32  -  1  )D 


-U)xUy((32  -  1) 
D2-n2-  u2y{f32  -  1) 

jUy(02  -  1  )D 


JUX(P2  -  1  )D 

jujy{(32  -  1)D 

02(D2  -n2)  +  Sl2f3‘ 


For  the  sake  of  simplicity  in  the  definitition  of  K  we  use  an  alternative  notation  for  strain: 


7 ij  =  2c  ij  for  i  zf. 
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with  fi2  =  w2  +  Uy  and  ;32  =  (A  +  2p)//t.  The  determinant  of  this  matrix  is 

det (KtCIQ  =  /3V  [D2  -  fi2]3 
and  this  system  therefore  has  non-trivial  solutions  of  the  form 

■Ui  =  ( Ai  +  BiZ  +  CiZ2)e~nz  i  =  x,  y,  z 

where  Ai,  Bi,  Ci  are  independent  of  2,  and  we  have  also  used  the  condition  that  displacements  stay 
finite  as  2  — >  00  to  include  only  terms  with  negative  exponents.  Substituting  this  back  into  (A. 3)  we 
obtain 

P2  +  V 


^  [  j U-’' x  A  X  T  jUlyAy  “h  BZ 


p2  - 1 


Bx  = 


'  n 


Bz 


By  ~  ~BZ 


which  gives 


Cx 

o' 

11 

=  CZ  = 

0 

Ux 

1 

0 

1 

JO  £ 

Ax 

Uy 

= 

e-nz 

0 

1 

jUJyZ 

n 

Ay 

ilz 

jwx 

.  n 

Q 

l/32  +  l  , 
n  p2  - 1  J 

.  S*  . 

(A. 4) 


The  three  still  unknown  terms,  Ax,  Ay,  B,  will  be  determined  by  boundary  conditions  at  2  =  0.  Here, 
we  will  first  pose  the  boundary  conditions  in  terms  of  the  load  distributions  on  the  surface  of  the  infinite 
halfspace,  defined  by  2  >  0.  The  2-normal  stress  azz  at  the  surface  2  =  0  equals  normal  load  on  the 
surface,  fz.  Similarly  the  shear  stresses  axz  and  ayz  at  2  =  0  equal  the  shear  loads  on  the  same  surface, 
fx  and  fy  respectively.  The  relevant  stresses  can  be  derived  by  substituting  (A. 4)  into  (A. 2)  which 
become 


Czz 


pe 


-f lz 


n 


“(2  +  a>2) 


jtoxCl 


2jujx 


T  H2 


J2-l 

~(w2  +  2 w2)  2 jojy  +  nz 

P2 


jiOyCl 


n 


p2  - 1 


-p  H2 


^1 

...  I 

Ay 

.  . 

(A.5) 


The  unknown  terms  Ax,  Ay ,  Bz  can  then  be  found  by  letting  2  =  0  and  either  invert  Eq.  (A. 4)  with 
[hi,  tty,  uz ]  =  tiyo,  u^o]  or  Eq.  (A.5)  with  [(Xx^,  &yzi  &zz\  =  [/11  fyi  /-]• 
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TFMs  for  Tractions  as  Inputs 


It  is  convenient  to  solve  for  each  load  case  separately  and  state  the  solution  in  terms  of  a  Transfer 
Function  Matrix  (TFM)  which  after  we  have  substituted  for  0  and  finally  becomes 


with 

p-zQ 

-  w 


2fl2  —  u>l(2v  +  flz) 
—ujxuiy(2u  +  flz) 

— ju>xf2(l  —  2v  +  flz) 


u  =  Tuff, 

—  U)XUy(2l'  +  flz) 

2fl 2  —  <jjy(2v  +  flz) 
—jujyQ(  1  —2  u  +  flz) 


ju>xfl(l  —  2v  —  Hz) 
jwy£l(l  —  2v  —  flz) 
(2(1  —  v)  +  flz)fl2 


(A.6) 


The  strain  and  the  stress  field  TFM’s  are  obtained  using  the  relations  in  (A. 2)  and  results  in 


with 


Trt  = 


2/tfl3 


with 


TV/  = 


fi3 


-  Teff,  where,  £  £yy,  ^zz>  cyz] 

ju>x  (2fl2  -  u>2(2i/  +  flz))  -ju^LOy  (2v  +  flz)  -flu;2  (1  —  2i/ -  flz) 

—jioxuSy  {2v  +  flz)  jiOy  (2fl2  -  u>y(2v  +  flz))  -flu;2  (1  -  2v  -  flz) 
-j?l2iox  (2v  -  flz)  -jfl2LOy  (2is  —  flz)  —  fl3  (1  —  2i;  +  flz) 

■jujy  (fl2  —  u>1(2v  +  flz))  juix  (fl2  -  u;2(2^  +  flz))  — flcjsuiy  (1  —  2v  —  flz) 
-fl2  (fl  —  w2z)  fl2wxu;j/z  jfl3LOxz 


fru^u;. 


2;U/y  ^ 


fl2wxu;j/z 
-fl2  (fl  -  urz) 


jn3u>y: 


(AT) 


and  <T  Tffff,  where,  (J  —  ^yy,  @zz,  0"xy,  &xzi  ^Tyz] 


ju>x  (2fl2  +  2^u;2  —  a2  zfl)  ju;y  (2z;u;2  -  u>2zfl)  -S1  (2wJ  +  w2(l  -  flz)) 


ju)x  ( 2vlo2x  -  flu;2z) 

jCl3LJxZ 


jcjy  (2fl2  +  2 vul  -  u>2zfl)  -fl  {2vu>l  +  u>2(l  -  flz)) 
jn3u)yz  -fl3  (1  +  zfl) 


juy  (fl2  -  u>l{2v  +  zf2))  jux  (fl2  -  u>2(2*;  +  zfl))  -flu;xu;y  (1  -  2v  ~  zfl) 


-fl2  (fl  -  u>2z) 

Q2U)xUlyZ 


tl2U)xUlyZ 

— fl2  (fl  —  u;2z) 


jQ3ivx: 


jfl3 


uaz 


(A.8) 
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We  also  calculate  the  TFM  for  the  mean  normal  stress  defined  as 


which  gives 


with 

r** 


Pn  =  -  ( crx  +(Ty+<Jz) 


Pn  — 


!Pi./ 


2 E  e~zCl  r 

T  2M7  L  juJx  juJy 


(A.9) 


The  relationship  between  uq  and  /  can  be  found  by  letting  z  =  0  in  Eq.  (A. 6)  obtaining 

uo  =  Tuoff, 


with 


Tuof  ~  2nW 


2(fl2  -  vtol) 

- <2l/UJXLJy 


-2vWxL0y 

2(Q2  -  VUJ2y) 


j(  1  -  2 v)<jJxSl 
j(l  —  2  v)u}y£l 


(A.10) 


[  —  j(l  —  2v)ujxQ  —j(l  —  2v)u}yCl  2(1  —  jy)02  J 

We  note  that  the  ;  component  decouples  from  the  x  and  y  components  when  the  material  is  incompress¬ 
ible  v  =  0.5. 


TFMs  for  Displacements  as  Inputs 


The  TFMs  when  the  surface  displacements  are  considered  the  inputs  are  found  similarily  as  the  ones  in 
previous  section: 

/  Tfuuu  o, 


with 


Tfua  (3  -  4 i/)fi 


SI2  (3  —  4^)  +  u)\ 

x  Wy 

2j(l  —  2is)u>xQ 


n2(3  -  Av)  +  u>‘ 
2j(l  —  2v)l Oyfl 


2j(l  —  2v)uxSl 
2j(l  -  2 v)u}yQ 

4(1  -j/)n2 


U  —  T'uut)  U0 , 


(A. 11) 
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with 


(3  —  4v)Sl 


(3  -  4i/)n  —  w2.z 

—  UJXUJyZ 

-ju!xSlZ 


- U!X(jJyZ 

(3  —  4 u)Sl  -  ujyZ 

—  jljJySlZ 


-jujxSlz 

—  jlVySlZ 

(3  —  4i/)Sl  + 


(A.12) 


For  the  strains 


~  —  FeU0TrQ, 


with 


e-n 

Tcu “  (3  -  4v)Sl 


M  ((3  —  4r y)Q-ujlz) 
-ju>xoj2yz 
—jSLox  (1  —  Slz) 

^iOy  ((3  -  4i')f2  -  2u>2z) 
|((3-4I,)n2+a>2(l-2nz)) 
—  ^U)XLOy(l  —  2zSl) 


-ju\UyZ 

jUy  ({2>  -  4v)Sl  -  u2yz) 

—  jSllxJy  (1  —  Slz) 

^UJX  ((3  -  4v)Sl  -  2 uj2z) 
-\ujxujy(l  -  2zSl) 
((3-4i/)ft2+w2(l-2fi  z)) 


flWyZ 

-SI2  (2(1  -  2v)  +  Slz) 

QuJXU)yZ 

jSlujx(l  -  2v)  +  Slz) 
jSlux(l  -  2v)  4-  Slz) 

(A. 13) 


and  for  the  stresses 

with 


T 

OXi 


wo, 


-  2  pe~^ 

au "  (3  -  4u)Sl 


x 


M  ((3  -  2u)Sl  -  u>2z) 
jwx  (2vSl  —  u j2z) 
—jSlwx(l  —  2v  —  Slz) 
i(jy  ((3  -  4v)Sl  -  2u2z) 

-  \  ((3  -  4v)Sl2  +  w2(l  -  2Slz)) 
—  \ujxLjy(  1  —  2Slz) 


jujy  (2 vSl  -w22) 
juiy  ((3  -  2v)Q  -  u2z ) 
—jSlu>y(  1  —  2v  —  Slz) 

\ux  ((3  -  4i/)ft  -  2 ui2yz) 
-±uxuy(l  -  2122) 

((3  —  4v)Sl2  +  w2(  1  -  217z)) 


Finally  for  the  mean  normal  stress  we  obtain 


— SI  (2vSl  —  lj lz) 

—SI  (2isSl  —  ix>2z) 

-SI2  (2(1  -  v)  +  Slz) 

Slid X  QJ y 

jfM(l  —  2v  +  Slz) 
jSlL0y(l  —  2v  +  Slz) 

(A. 14) 


Pn  —  TPnUoUo, 
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with. 


L  p™  ^0 


2  E  e~zn  r 


n 


(A.15) 
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Appendix  B 


Effects  of  Finite  Sensor  Size  on  a 
Stress  Field 


If  the  sensor  is  of  a  finite  size,  it  will  disturb  the  stress  field  around  it  and  one  may  need  to  corrects  for 
these  effects  at  that  sensor.  There  is  also  the  possibility  that  the  sensor  will  disturb  the  sensor  field  in 
such  a  way  that  it  will  effect  other  sensors  close  by.  It  is  therefore  appropriate  to  assess  how  close  sensors 
of  finite  size  can  be  placed  without  disturbing  each  other.  We  do  that  by  investigating  how  a  spherical 
inclusion  with  different  elastic  properties  than  the  medium  around  it  reacts  to  an  external  stressfield. 


We  assume  that  the  stress  field  is  uniform  away  from  the  inclusion  or  so  slowly  varying  that  it  can 
be  assumed  to  be  uniform.  It  then  follows  from  dimensional  analysis  that,  under  classical  elasticity 
assumptions,  the  stress  intensity  factor1  is  independent  of  the  absolute  size  of  the  inclusion  and  only 
depends  on  its  geometrical  form  [24]. 


Here  below  we  give  the  solutions  for  a  cavity2  under  3  different  loading  conditions,  i.e.  under  (i) 
uniform  stress,  (ii)  uniaxial  stress,  and  (iii)  pure  shear. 


1  The  stress  intensity  factor  is  defined  as  the  ratio  of  the  maximum  stress  to  the  uniform  stress  at  a  distance. 

2If  t  he  inclusion  is  not  a  cavity  but  an  inclusion  of  a  different  material,  the  behavior  will  be  quite  similar  with 
slightly  lower  stress  intensity  factors  [24]. 


141 


Uniform  Stress 


Figure  B-l:  Shows  the  effects  of  an  inclusion  on  a  uniform  stress  field.  The  plot  shows  the  radial  stress 
given  by,  aRR  =  (l  -  a3/R3). 

1.  Uniform  Stress 


In  the  case  of  a  uniform  stress,  the  radial  stress,  aRR  and  the  hoop  stresses  erg#,  are 

-3  \  /I  a3 


&RR 


1 


a' 

R3 


and 


—  &06  —  CTqo  |  1  “h  gy  j 


which  means  a  stress  intensity  factor  of  3/2  for  the  hoop  stresses.  The  radial  stress  is  plotted  in  Fig. 
B-l  and  there  we  see  that  the  error  is  4%  at  R  —  3 a,  and  1%  at  R  =  5 a  where  R  is  the  distance  form  the 
center  of  the  inclusion  and  a  is  its  radius.  By  that  measure,  sensors  would  have  to  be  put  5  diameters 
apart  for  less  than  2%  cross-sensor  error. 


2.  Uniaxial  Stress 


In  the  case  when  the  external  loading  is  a  uniaxial  tension  a <*,  applied  at  infinity  in  the  ^-direction 
then  the  most  significant  variations  are  observed  in  the  normal  stress  on  the  plane  perpendicular  to  the 
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Uniaxial  Tension 


2 


Figure  B-2:  Shows  the  effects  of  an  inclusion  on  an  uniaxial  stress  field.  The  plot  shows  the  --normal 
stress  given  by,  aZ2  =  (l  +  a3 /(6R3)  +  a5/R5). 

loading,  i.e.  azz,  which  becomes 

/  14-5 v  a3  9  1  a5  \ 

°zz  “  a°°  V  +  2  7  -  5u  R3  +  2  7  -  5 v  R5  ) 

=  +  H  +  to  ^o-5 

which  means  a  stress  intensity  factor  of  13/6.  The  azz  stress  is  plotted  in  Fig.  B-2  and  there  it  can  be 

seen  that  the  error  is  1%  for  R  =  3 a  and  0.2%  for  R  =  5a.  The  effect  on  the  stressfield  are  therefore  not 

spreading  as  far  as  in  the  case  of  a  uniform  stress. 


3.  Pure  Shear 


If  uniform  tension  is  applied  in  the  ^-direction  and  and  uniform  pressure  is  applied  in  ^-direction  the 
net  effect  is  pure  shear  in  the  planes  x  =  ±z,  ie.  axz  =  a0 0  [49].  It  has  a  maximum  at  x  =  z  =  0,  y  =  ±a 
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Pure  Shear 


z 


HHIIUJI 


Figure  B-3:  Shows  the  effects  of  an  inclusion  on  pure  shear.  The  plot  shows  the  arz-shear  stress  along 
the  y- axis,  given  by  axz  —  (l  +  2a5/(3i?5)).  We  observe  that  the  disturbance  is  more  localized  than 
the  ones  in  Figs.  B-2  and  B-l 


and  along  the  y- axis  the  solution  is 

az.  = 


which  means  a  stress  intensity  factor  of  3/2.  The  axz  stress  is  plotted  in  Fig.  B-3  and  there  it  can  be 
seen  that  the  error  is  only  0.3%  for  R  =  3 a  and  0.02%  for  R  —  5 a.  These  results  indicate  that  under  the 
assumptions  given,  shear  sensors  are  less  sensitive  than  normal  stress  sensors  to  disturbances  caused  by 
their  own  finite  size.  We  note  that  this  observation  is  analogous  with  the  observations  made  in  section 
1  regarding  the  lowpass  properties  of  the  mean  normal  stress  (spread  out)  and  the  bandpass  properties 
of  the  shear  strain  (localized). 


5(1  -  2 z/)  a3 
1  -  -4 — r-^—  + 


a'-' 


1  + 


7  —  oi/  R3 

2  a 5_ 

3  R5 


7  —  5i/  R5 
for  v  —  0.5 


In  summary  the  results  of  the  above  are  that  the  disturbances  of  the  stress  field  caused  by  finite 
size  sensors  estimated  using  solutions  of  spherical  inclusions  in  an  infinite  medium  have  the  following 
characteristics: 


•  The  magnitude  of  the  effect  depends  only  on  the  geometry  of  the  sensor,  not  on  its  dimension. 

•  Sensors  from  a  material  that  is  stiffer  than  the  surrounding  medium  have  generally  a  lesser  effect 
than  sensor  made  from  softer  materials. 
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•  The  effects  on  the  mean  normal  pressure  die  out  oc  ( a/R )3  where  a  is  the  radius  of  the  sensor  and 
R  is  the  distance  from  its  center. 

•  The  effects  on  shear  stress  die  out  as  oc  (a/R)5,  and  shear  sensors  are  therefore  less  likely  to  be 
sensitive  to  the  effects  of  nearby  sensors. 
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Appendix  C 


Derivation  of  Dynamic  Transfer 
Function  Matrices 


In  order  to  solve  the  nine  partial  differential  equations  given  by  Eqs.  (4.1)  and  (4.2)  in  chapter  4  section 
2.1,  we  make  use  of  the  triple  Fourier-Laplace  transform  defined  by 


/**fOO  p  p  +  OO 

f{ujx,ujy,s)=  //  f{x,y,Ti)e~3(u'x+u>''y)+STl&x&y&Ti 

Jo  J  J  —  oc 

i  r-f+joo  r  p+ oo 

5,2/,Tl)=8 J  jj  f(ux,uy,s)eAu’rX+“'’y)+STldLoxdujyds. 


(C.l) 


Assuming  that  all  the  integrals  are  well  behaved  as  \x,y,T\\  — *  oo,  we  apply  this  transform  to  each 
term  in  the  equations  in  (4.1)  and  (4.2).  Then  the  operator  ^  is  replaced  by  multiplication  by  jujx, 
similarly  ^  is  replaced  by  multiplication  by  ju)y.  We  further  use  D  to  notate  jp  and  assuming  zero 
initial  conditions  at  t  =  0,  we  replace  ^  by  s.  The  9  partial  differential  equations  have  now  been 
transformed  into  ordinary  differential  equations  in  D.  ( W  becomes  ^).  We  write  (4.1)  and  (4.2)  using 
matrix-vector  notation  as 

I<Td  =  02s2u 


and 


where 


Kl 


jux 

0 

0 


a  =  Ce  e  =  Ku 
0  0  ju)y 

jujy  0  jujx 


D 

0 


0 

D 


0  D  0  jU>X  jiOy 
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c  = 


A  -t-  2/x  A  A  0  0  0 

A  A  4-  2  fi  A  0  0  0 

A  A  A  -f-  2  fz  0  0  0 

0  0  0  fL  0  0 

0  0  0  0  /X  0 

0  0  0  0  0  /x 


O  ®yy>  &zz>  & x\ /i  &xzi  &yz]  ?  ti  —  [^x?  t£y,  ^ z ]  ■>  &tld  £  —  [^xu  £ yyi  &zzi  'Ixyt  ~fxz)  "7 yz 

nating  a  and  £  we  get 

\KtCK  -  02s2l]  u  =  0 


*  ir  i 


Elimi- 


(C.2) 


where 


K1  CI< 


D2  -  O2  -  w2(/?2  -  1)  -u>xuy((32  -  1)  jux(02  -  1  )D 

-UXLOy(02  -  1)  D2  -Q.2  -  U>1((32  -  1)  jLOyiP2  -  1  )D 
jwx(02-l)D  My(P2-l  )D  p2(D2-n2)  +  n202 


with  fi2  =  u>l  4-  u!2  and  02  =  (A  +  2 n)/ii.  The  determinant  of  this  matrix  is 


det {KtCI(  -  ft2 s2 1)  =  /?V  ( D 2  -  712))  (D2  -  n22) 


where 

nx  =  (fl2  +  s2)1/2  ri2  =  (fl2  +  /32s2)1/2 
This  system  therefore  has  non-trivial  solutions  of  the  form 

m  =  Aie~niz  +  (Bi  +  Ciz)e~nsz  i  =  x,  y,  z 


where  we  have  used  the  condition  that  displacements  stay  finite  as  z  — >  oo  to  include  only  terms  with 
negative  exponents.  Substituting  this  back  into  (C.2)  we  further  get  that 


A. 

=  ±(jv,A, 

Bx 

juJx 

~  sTs= 

By 

ju>y  n 

-  -~nB‘ 

cx 

II 

£ 

II 

P 

1 1 

P2  +  l\ 
02-l) 


1 


For  the  sake  of  symmetry  we  use  the  equivalent  definition  of  strain:  7 ij  =  2s ij , 


*  ^  j- 
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which  gives 


ux  =  -AzJ-^e~n'z  +  Bxe~n*z 

Til 

Uy  =  -Azi^-e~n'z  +  Bye-n2Z  (C.3) 

7l\ 

uz  =  Aze~niz  +  —  ( u>xBx  +  joJyBy)  e~n2Z. 


The  three  still  unknown  constants,  Az,  Bx,  By  will  be  determined  by  boundary  conditions  at  2  =  0. 


TFMs  with  Tractions  as  Inputs 


We  will  first  pose  the  boundary  conditions  in  terms  of  the  load  distributions  on  the  surface  of  of  the 
solid,  defined  as  the  infinite  halfspace,  defined  as  2  >  0.  The  2-normal  stress  azz  at  2  =  0  equals  normal 
load  on  the  surface,  2  =  0,  fzo  and  similarly  the  shear  stresses  axz  and  ayz  at  2  =  0  equal  the  shear 
loads  on  the  same  surface,  fXo  and  f,JU  respectively.  The  relevant  stresses  can  be  derived  by  substituting 
(C.3)  into  (4.2)  and  then  become 

n2 

A~2jujye-n'z  -  UxUJyB*  +  (u2.  +  u^Bv  e-n»* 

n-2 

2f)2  4.  s2/?2 

-Az - ±±JLe-m:  _  2 j(LoxBx+uyBy)e~n2Z. 

n  1 

It  is  convenient  to  solve  for  each  load  case  separately  and  state  the  solution  in  terms  of  a  Transfer 
Function  Matrix  (TFM)  which  finally  becomes 


Gxz 

TTyz 

M 

<Tzz 

t1 


uj2xh2 

juxn\ 

Tfu  = 

U)XWyTl2 

u2yn  2 

ju>ynj 

e~niz 

2hF2 

_  jujxnin2 

juynin2 

-nml  _ 

ujy{2nin2  —  n|)  —  nlnl 

ojxu>y(n^  —  2nin2) 

jujxrnn\ 

+ 

uxwy(n%  -  2nin2) 

u>l(2niii2  —  n\)  —  n^n22 

ju>yniTi2 

e~n2Z 

2  iin2F3 

-jujxn2nl 

-jujyn2nl 

rii  n2Q.2 
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where  n3  —  (fi2  +  /?2/2 s2)1/2  and  F3  =  n\  -iiin2Q2.  The  strains  are  obtained  using  Eq.  (4.3)  and  hence 
the  TFM  is 


-ni  J 


2  HF3 


juln2 


jLoxu>‘n2 

-juxn\n2 

j2uj2u)yn2 

-2w2nin2 


LO^yn2 

ju3yn2 

-juyn\n2 

j2ojxbJ2yn2 
2ujx^'y  Ti\Ti2 


2„2 


-2ujxLOynin2  -2a;2n1n2 


z'l3 

n\nl 

—  2tl)  XOJyTl2 

-j2uxrnnl 

-j2ujyninl 


+ 


2  fm2F3 


ju>x(2u)2nin2  -  n2(2n2  -  w2)) 

jtJluy(n\  -  2n,\n2) 

ju>x'J2(n'l  -  2n1n2) 

jtoy(2uln1n2  -  n§(2n§  -  u>2)) 

w27iin| 

ju>xnln\ 

j{j}yn\n$ 

-nin2Q2 

j 2wj,(ni ?i2(w2  -  w2)  +  n|(w2  -  n|)) 

j2ux(nln2{J2x  -  w2)  +  nf  (w2  -  n§)) 

2wxWj,nin.2 

2  n2  n|  —  2cj^ni  n  2 

2uxujynin\ 

j2u>xn\n2n\ 

2n2n\  —  2w2n!Ti2 

j2uynin2n2 

(C.5) 


and  the  stresses  are  obtained  using  Hooke’s  law  in  Eq.  (4.2),  hence 


Tfa  = 


2  F3 


-ju>xn2((P2  -  2 )s2 

-2u£) 

-juyn2({P2  -  2)s2  -  2u>2x) 

~((P2  -  2)s2  -  2 Jinf) 

-juxn2{{P2  -  2)s2 

-H) 

-ju)yn2({(32  -  2)s2  -  2n>2) 

~((P2  -2)s2  -2uj2ynl) 

-2ju>xn2n2 

—2juyn2n2 

2  n3 

j2J1xuyn2 

j2uxw2yn2 

-2  WjWy/ll 

-2^71x712 

—2u>xLOynin2 

—j2ujxnin^ 

—  2  71x712 

-2u2Tiin2 

—j2ujynin2 

+  2  n2F3X 

j2iox{2uj2ynln2  -  ra§(2n§  -  u>2)) 

-  2?iin2) 

2  jujxri22nl 

j2wy(n1n2(ul  -  w2)  +  nf(w2  -  n§)) 
2n2  Jig  —  2w^niTi2 

2u)xuJynin2 


j2LU2xu)y(n23  -  2njn2) 

2w2nin| 

j2u)y(2uj'2n1n2  -  n|(2n|  -  w2)) 

2w2n1?i2 

2 juynlnl 

—2niri2^l2 

j2u>x{n1n2(uj2x  -  w2)  +  n§(a;2  -  n§)) 

2u)xu}ynin2 

2u>xujyninl 

j2ioxri\n2n2 

2n2n\  —  2(jj2n\n2 

j2ujynin2n\ 

The  mean  normal  stress  defined  as  pn  =  (axx  +  ayy  +  <rzz)/3  is  then 


/p.. 


-jLOxii-i  -jujyn2  n2 


3 P2  —  4  s2e  ni‘ 

3 


(C.6) 

(C.7) 


TFMs  with  Displacements  as  Inputs 


When  the  displacement  distribution  on  the  surface  is  taken  as  inputs,  we  use  Eq.  (C.3)  to  solve  for 
(.4,,  Bx,  By)  which  gives: 


The  TFM  giving  the  subsurface  displacements,  u  =  TUU(]uo,  as  being 


r  w2 

X 

CDxa;y 

jwxn2 

g  Tli  2  _  e~n22 

‘  1 

0 

0  ' 

U)XUy 

Id2 

y 

ju>yn2 

fi2  -  «in2 

0 

1 

0 

jUJyTll 

nin2 

_  0 

0 

1  _ 

(C.8) 
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the  TFM  for  the  strains,  e  =  Tue uq  ,  as  being 


+ 


JW2  j^Uy 

-w2n2 

ju)xL>l  ju^l 

-w2n2 

-juxn\  -jvyn2 

n2n.2 

e 

•n\z 

2juj2xWy  2juxu)l 

-2ujxwyn2 

n2  - 

■  nin2 

— 2uilni  —2u)xu>yni  — 2jcnxni7i2 

2ojx^yii\  2uiyn  x 

-2juymn2  j 

jux( L02y  -  nin2) 

-j^luy 

u2xn2 

-juxul 

juy(uj2  -  n1n2) 

w2n2 

juxnxn2 

jujynin2 

02b2 

M(w2  -  w2  -  nin2) 

jwx(w2  -  w2  -ni?i2) 

2u1XU)y1\2 

ni(n?2  —  io2)  —  ui2n2 

wxwy(ni  +  n2) 

jojx{  ft2  +  n2) 

uxLjy(n  i  +  n2) 

nx{n\  -  w2)  -  w2n2 

;u;y(n2  +  n|) 

(C.9) 


n2  —  n in  2 ’ 


the  TFM  for  the  stress  a  =  Tua  uq  ,  as  being 

jojx(  2u/2  +  s2((i2  +  2))  jioy(2ujl  +  s2((32  +  2)) 
j'u>x(  2ujy  +  s2((32  +  2))  jujy(2uj2  +  s2(i32  +  2)) 


T  — 

-'■0'uo  — 


-n2(2Sx+s2(P2+2)) 
-n2(2uj2y+s2(P2  +  2)) 


+ 


2jWx^3 

ZjUyU'i 

-2 n2n2 

fte~ni: 

2j^xu)y 

2  juxw2y 

—2  ujxu>yn2 

u2  —  ??  i  ?r2 

~2u)2xm 

-2ju>xnin2 

-2  LOxLOylll 

-2w2ni 

-2ju>yiiin2 

2 i^rr(w2  -  nin2) 

-2jJluy 

2w2n2 

-2jwx.u2 

2ju)y{ul  -nin2) 

2w2n2 

2  jojxiiin2 

2jojyniii2 

2  fi2n2 

pe~n2Z 

jux{u>l  -u2x-  iiin2) 

juy(ul  -o>2  -nin2) 

2'jjxu)yn2 

fi2 —  nin2 ’ 

2nJ?ii  -  w2(rii  +  n2) 

uxu;y(n  i  +  n2) 

jujx(Q2  +  n\) 

2w2tii  —  w2(?ii  4-  n2) 

juy(Cl2  +  nl) 

and  finally,  the  TFM  for  the  mean  normal  stress,  pn  —  TPnUuu0,  as  being 

lie~n'z  3/32  —  4 


T  = 

X  P„  K(J 


JUJX  Jljjy  n-2 


fi2  —  n\n2  3 


(C.10) 


(C.ll) 
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The  variable  s  in  the  above  equations  is  the  Laplace  variable  corresponding  to  the  scaled  t\  time 
variable.  To  get  the  solution  in  terms  of  a  Laplace  variable  corresponding  to  the  actual  time,  t,  we  need 
only  to  replace  .s'2  everywhere  with  c~2s'2  so  that  we  will  get  (after  we  drop  the  prime) 

ni  =  \/n2  +  Cp2s 2  ns  =  \/f!2  +  cj2s2  n3  =  Jtt2  +  Cp2s2/ 2. 
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Appendix  D 


Retrival  of  the  Static  Solution 


The  dynamic  solution  should  approach  the  static  solution  when  s  — +  0,  but  since  both  the  denominator 
and  the  numerator  vanish  for  s  =  0  we  use  Taylor  expansion  for  the  denominator  and  also  the  (e~n2Z  — 
e~UlZ)  term  around  s  =  0.  Writing  each  step  out  explicitly  we  get 

ni  =  \Jfl2  +  Cp2s 2  =  0(1  +  |cp2s2)  +  0(s4) 
n-2  =  \/02  +  cj2s2  =  0(1  +  ±cj2s2)  +  0(s 4) 


and  hence  the  denominator  becomes 


F3  =  (O2  +  |cs/l2s2)2  -  027ii?i2 

=  O4  +  0 2c~2s2  -  02(02  4-  \c~2s2  +  \c~2s2)  +  0(s4) 
=  hn2s2(cl 2  -  c~2)  +  0{s4), 

and  the  numerator  term 


(, e~n2Z  -  e~n'z) 


=  exp  (-Qz(l  +  ±cs2s2n  2))  -  exp  (-Oz(l  +  ^cp2s2Cl  2))  4-  0(s4) 
—  e~Qz  [exp  (-|c7ft2s20-12)  -  exp  (-|cp2s20-1z)]  4-0(s4) 

=  e“n*  (1  -  |c7h2s20-12  -  (1  -  Icp-VO-1*))  +  0(s4) 

=  |e-n2s20-12(c-2  -cjh2)  +  0(s4). 
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Substituting  this  into  Eq.  4.5,  we  get  for  the  mean  normal  stress 


and  the  shear  strain 


lims — ,0 


3/32  -  4  s2  (A2  +  hc7hs2)e~n'Z 
W2  clh  fi2s2(c7h2  -  Cp2) 

iff  ~_4  e-n« 

3(/?2-l) 


lim5_0  re„/. 


1  ja;xn3  (§s2Q-1e-f2z^(c~2  -cM2)) 

|n2s2(c;h2-Cp-2) 
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Appendix  E 


The  Kalman  Filter  With  Static 
Load 


Since  the  Wiener  filter  given  in  chapter  3  and  the  Kalman  filter  given  in  chapter  4,  minimize  the  same 
error  measure,  they  should  become  identical  when  the  (loads)  signals  become  static.  This  can  be  seen, 
most  easily,  by  using  the  discontinuous-time  version  of  the  Kalman  filter  under  the  assumption  that  the 
input  is  an  unknown  constant.  In  this  case  the  Kalman  filter  estimate  is  given  by 


£  = 

I<  = 

PC*R~1 

-l  _ 

Po'1  +C*R~lC 

(E.l) 


We  identify  C  as  the  static  TFM  T(uix,  wy,  z),  since  it  is  the  encoding  solution  relating  the  input  and  the 
state  vector,  which  only  contains  the  sensor  output  under  static  conditions.  Similarily  we  note  that  in 
the  static  case  the  sensor  noise  intensity  matrix  R(ojx,u>y)  is  the  power  spectrum  matrix  for  the  sensor 
noise,  Vn(wx, uy),  also  the  initial  covariance  estimate,  P0(uXlu)y)  is  the  power  spectrum  matrix  for  the 
solution  Ps{ujx,u)y),  or  R  =  Vn  and  Po  =  Vs-  Substituting  the  above  into  the  equations  in  Eq.  (E.l) 
and  combining  we  finally  obtain 

|  =  [Ps1 +T*P^1T]-1T'V^1  C 
=  VS[VN  +T*VsT}~lT*  C 
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which  we  recognize  as  the  Wiener  filter  derived  in  [38,  80]  and  used  for  the  static  case  in  chapter  3. 
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Appendix  F 

General  Statement  of  the  Decoding 
Problem 


Problem  Statement:  Given  a  3D  sensor  signal  5  =  [SI,  52,  53]T,  a  function  of  the  strain  tensor  at 
the  sensor  location,  what  conditions  does  it  have  to  satisfy  in  order  to  be  suitable1  for  decoding. 


Remark:  This  approach,  of  assuming  a  3D  loading  and  a  3D  sensor  signal  makes  it  necessary  to  consider 
strain  measures  in  sets  of  three  and  one  can  therefore  not  answer  questions  such  as:  Is  strain  energy 
density  suitable  as  relevant  stimulus.  (One  will  need  to  specify  two  additional  strain  measures.) 


Assumptions 


A1  The  sensors  form  a  continuous  2D  sheet  that  is  located  in  some  sense  parallel  to  the  contact 
surface  and  both  are  anchored  at  an  outer  boundary  where  boundary  conditions  are  known.  (The 
boundary  is  at  infinity  for  the  infinite  halfspace.) 

A2  We  further  assume  that  the  displacements  and  strains  are  parametrized  everywhere  by  three 
variables  ( x,y,z ),  z  being  constant  for  both  the  surface  and  the  sensor  sheet  (depth).  We  can 
then  write  the  sensor  signal  as  5  =  S(ui,uitj)  where  i,j  =  x,y,  z  and  u^j  is  the  partial  derivative 
of  xn  wit.  j  (e.g.  dux/dz  =  rtx>2). 


1  By  'suitable  for  decoding’  we  mean  that  there  is  a  one-to-one  mapping  between  all  allowable  traction  distri¬ 
butions  and  sensor  distributions. 
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A3  Finally,  we  assume  that  the  sensor  sheet  is  sufficiently  close  to  the  surface  so  that  there  is  a 
one-to-one  relationship  between  surface  displacements  and  the  ones  at  the  sensor  depth. 


Conditions:  A  sensor  signal  will  be  suitable  for  decoding  if  it  satisfies  the  following  conditions: 


Cl  The  relation  S  =  S(ui,Uij)  i,j  =  x,y,z  can  be  integrated  wrt.  z  using  the  encoding  solution2 
so  that  it  becomes  5' =  5'(ui,Uitfc)  i  =  x,y,z  k  =  x,y. 

C2  The  relation  S'  =  S'(ui,Uitk)  is  integrable  wrt.  x,  y  giving  S"  =  S"(ui)  i  =  x,y,z. 

C3  This  relation,  S",  has  to  invertable,  i.e.  there  exists  a  function,  5"_1 ,  such  that  u  =  [ur ,  uy,  uz]T  = 


Remark:  Condition  C2  is  different  from  Cl  as  does  not  need  the  encoding  solution,  but  utilizes  the 
fact  that  the  integrand  is  known  over  all  the  integration  domain  and  at  the  boundary.  It  is  therefore 
a  condition  on  the  sensor  signal  independent  of  the  particular  shape  or  mechanical  composition  of  the 
fingerpad. 


2By  ’’the  encoding  solution”  we  mean  the  model  of  the  fingerpad  that  given  the  displacements  at  one  depth 
gives  the  displacements  and  also  strains  at  all  depths  and  hence  the  sensor  signal  S  for  all  z. 
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